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Abstract 

Struetural  health  monitoring  (SHM)  is  a  eompetitive  approaeh  for  damage  deteetion  in  aireraft 
struetures,  wherein  online  information  is  eolleeted  and  eompared  with  an  existing  database  for 
the  undamaged  strueture,  to  obtain  real-time  information  about  the  presenee  of  damage.  The  goal 
of  this  researeh  is  to  develop  numerieal  models  of  inverse  problems  for  damage  deteetion  in 
aireraft  struetures,  whieh  eould  later  be  part  of  an  on-board  system  for  SHM.  In  this  work,  the 
numerieal  modeling  has  two  main  branehes:  I.  The  direet  problem:  a  model  is  required  to  obtain 
information  on  the  distribution  of  the  quantity  of  interest  throughout  a  given  damaged  strueture. 
The  model  of  the  direet  problem,  using  the  boundary  element  method  (BEM),  is  expeeted  to 
reproduee  the  reality  of  an  aireraft  strueture.  II.  The  inverse  problem:  a  model  is  required  to 
loeate  the  struetural  damage  given  the  information  on  the  quantity  of  interest  at  partieular 
loeations  (sensor  loeations).  To  inerease  the  reliability  of  the  deteetion  approaeh,  a  eombination 
of  independent  optimization  and  identifieation  proeedures  ean  be  used.  Some  treatment  of  the 
model  uneertainties  is  required,  due  to  the  stoehastieity  in  the  problem  variables  and  parameters. 


1.  Introduction  /  research  outline 

1.1  Background 

Aireraft  struetures  are  subjeet  to  damage  during  their  useful  life.  The  timely  deteetion  of 
damages  in  aireraft  struetures  is  an  important  feature  for  flight  safety.  The  usual  inspeetion 
proeedures  during  regular  maintenanee  intervals  may  lead  to  problems  sueh  as: 

•  Inspeetion  intervals  might  be  too  large,  thus  allowing  damage  to  propagate  unnotieed  for  an 
unaeeeptable  time  interval  or  through  an  unaeeeptable  extension; 

•  Critieal  struetural  eomponents  might  be  diffieult  to  aeeess,  thus  imposing  disassembly  / 
assembly  proeedures  whieh  are  time-eonsuming  and  expensive,  sometimes  requiring  jigs  and 
other  speeial  tools  for  proper  assembly  of  the  strueture; 

•  Some  non-destruetive  teehniques  (sueh  as  eddy  eurrent,  for  example)  may  be  portable,  not 
requiring  full  disassembly  of  the  strueture,  but  might  be  inaeeurate  or  might  depend  strongly 
on  the  teehnieian’s  experienee  for  a  damage  to  be  deteeted  properly; 

•  Other  non-destruetive  teehniques  (sueh  as  X-Ray,  or  magnetie  partieles,  for  example)  might 
require  full  eomponent  disassembly  and  removal  to  an  industrial  faeility  to  perform  the 
struetural  inspeetion. 

•  Some  struetural  eomponents,  for  example  those  made  of  eomposite  materials,  may  present 
internal  damage  whieh  is  diffieult  to  deteet  using  standard  inspeetion  teehniques,  as 
ultrasonie  tests. 

1.2  Scientific  Challenge 

Struetural  health  monitoring  (SHM)  is  a  eompetitive  approaeh  for  damage  deteetion,  wherein 
online  information  is  eolleeted,  eompared  with  an  existing  database  for  the  undamaged  strueture, 
and  from  this  eomparison,  real-time  information  about  the  presenee  of  damage  is  obtained,  its 
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loeation,  length,  speed  of  propagation,  and,  ultimately,  the  remaining  operational  life  of  the 
struetural  eomponent. 

Some  ehallenges  related  to  an  effieient  on-board  struetural  health  monitoring  system  inelude: 

•  The  system  must  be  small  in  size  and  weight,  must  eonsume  a  small  amount  of  power,  and 
must  not  interfere  with  the  aireraft  eleetrieal  system; 

•  The  system  must  be  reliable,  and  the  information  on  the  loeated  damage  must  also  be 
reliable:  thus,  the  system  must  have  redundaneies  in  the  built-in  numerieal  eodes,  whieh  must 
be  based  on  separate  independent  numerieal  models. 

1.3  Objective 

To  study  and  develop  numerieal  models  of  inverse  problems  for  damage  deteetion  in  aireraft 
struetures,  whieh  eould  later  be  part  of  an  on-board  system  for  struetural  health  monitoring. 

1.4  Approach 

The  approaeh  in  this  researeh  work  is  to  investigate  the  feasibility  of  numerieal  proeedures  for 
damage  deteetion  in  aireraft  struetures.  The  numerieal  modeling  eonsists  in  two  main  parts: 

I.  The  direct  problem:  a  model  for  the  strueture  is  required  to  obtain  information  on  the 
distribution  of  the  quantity  of  interest  (for  example,  the  aeoustie  pressure  or  the  stress  field) 
throughout  the  strueture,  given  the  boundary  eonditions  and  the  presenee  of  the  damage. 
The  modeling  of  the  struetures  is  earried  out  using  the  boundary  element  method  (BEM). 
The  advantages  of  the  Boundary  Element  Method  over  other  numerieal  methods  (sueh  as 
finite  elements)  are  well  known  from  the  literature,  espeeially  for  the  treatment  of  high 
gradient  problems,  sueh  as  the  stress  gradient  due  to  eraeks.  Numerieal  models  for 
potential,  aeousties,  or  elastieity  ean  be  used  in  eombination  or  independently,  to  simulate 
the  multiple  physies  present  in  lamb  waves,  stress  waves,  aeoustie  emission,  ete,  involved 
in  the  usual  struetural  monitoring  teehniques  ([1]  -  [3]). 

II.  The  inverse  problem:  a  model  is  required  for  the  proeedure  of  loeating  the  damage  in  the 
strueture  given  some  (partial)  information  on  the  quantity  of  interest  (for  example,  the 
stress  field)  at  some  partieular  loeations  (for  example,  where  some  sensors  are  plaeed).  Eor 
this  inverse  problem,  both  optimization  proeedures  (loeal  and  global  optimization)  and 
identifieation  teehniques  ean  be  used. 

Also,  both  the  direet  and  the  inverse  problems  are  in  faet  stoehastie,  and  involve  some  level  of 
treatment  of  the  randomness  of  the  parameters  and  variables  of  the  models. 

The  researeh  projeet  eoneentrates  on  three  main  problems: 

1.  The  direct  problem:  the  model  of  a  strueture  with  a  known,  assumed  damage.  Eor  this 
problem,  the  numerieal  eodes  to  investigate  must  inelude  the  model  of  reinforeed  panel 
struetures,  both  metallie  and  eomposite.  Some  possibilities  for  these  models  are: 

•  The  elastie  modeling  of  eraeked  anisotropie  plates  using  boundary  element  methods. 
Eurther  study  might  inelude  elastoplastie  and  elastodynamie  behaviors; 

•  The  aeoustie  modeling  of  damaged  anisotropie  plates  using  boundary  element 
methods.  The  study  might  inelude  aeoustie  propagation  from  a  generated  signal  or 
from  existing  aerodynamie  noise. 
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•  The  dynamie  modeling  of  the  vibration  signature  of  the  damaged  strueture,  ineluding 
natural  frequeneies  and  modes  of  vibration. 

2.  The  inverse  problem:  a  global  optimization  proeedure  to  minimize  a  funetional,  obtained 
from  differenees  between  measured  values  and  values  generated  by  the  numerieal  eode  at 
different  assumed  damage  loeations.  The  minimum  value  of  the  funetional  will  oeeur  when 
the  distanee  between  the  real  damage  and  the  assumed  damage  loeation  in  the  numerieal  eode 
is  also  a  minimum,  thus  giving  an  indieation  of  the  damage  loeation  and  size.  Some 
possibilities  to  model  this  optimization  problem  inelude  heuristies  sueh  as  evolutionary 
algorithms  (gene tie  algorithms,  or  differential  evolution,  for  example).  Additionally  to  the 
optimization  proeedures,  identifieation  teehniques,  sueh  as  artifieial  neural  networks  (ANN), 
ean  also  be  used  for  this  inverse  problem,  by  setting  the  desired  loeation  and  size  of  the 
damage  as  parameters. 

3.  The  modeling  of  uncertainties:  some  variables  involved  in  the  proeess  (sueh  as  the 
aerodynamie  loads),  do  not  present  determinist  values  and  must  be  treated  as  random 
variables.  Also,  some  parameters  of  the  strueture,  sueh  as  the  elastie  properties  and 
eonstitutive  behavior,  are  also  non-deterministie  and  must  be  identified.  The  treatment  of  the 
stoehastie  nature  of  the  problem  leads  to  parameter  identifieation  proeedures  (sueh  as 
Kalman  filter  identifieation,  for  example)  and  to  stoehastie  optimization  proeedures  (sueh  as 
response  surfaee  methodology  or  Monte  Carlo  simulation).  Proeedures  to  obtain  the  response 
surfaee  might  inelude  design  of  experiments  eombined  with  regression,  or  the  learning  of  the 
struetural  behavior  through  a  neural  network  proeedure. 

The  researeh  eovered  a  three-year  period  (from  2006  to  2009),  as  a  eollaborative  effort  between 
researehers  working  in  the  eomputational  meehanies  area  from  UNIFEI  (Federal  University  of 
Itajuba)  and  UNICAMP  (State  University  of  Campinas),  both  in  Brazil.  During  the  three  years  of 
the  projeet,  the  researeh  also  lead  to  advising  graduate  students,  both  at  UNICAMP  (with  the 
researeh  eoneentrated  in  the  study  of  the  behavior  of  eraeked  eomposite  plates  using  boundary 
element  methods)  and  at  UNIFEI  (with  the  researeh  eoneentrated  in  the  use  of  deterministie  and 
stoehastie  optimization  and  identifieation  teehniques  for  a  given  direet  modeling  (elastieity  and 
aeousties,  for  example). 

The  modeling  of  uneertainties  was  ineipient  in  this  work,  and  is  an  on-going  effort,  whieh  is 
planned  to  eontinue  as  a  eollaborative  researeh  work.  Also,  a  eombination  of  EEM  and  BEM  as 
direet  models  was  not  eovered  in  this  researeh  projeet,  and  is  being  planned  as  an  on-going 
researeh  work,  too. 

1.5  Resources  /  Research  team 

•  Professors  Ariosto  Bretanha  Jorge  (PI  -  Prineipal  Investigator)  and  Sebastian  Simoes  da 
Cunha  Jr.  (UNIEEI); 

•  Professors  Paulo  Sollero  and  Eder  Eima  de  Albuquerque  (UNICAMP); 

•  Students  from  the  Computational  Meehanies  groups  (both  from  UNIEEI  and  UNICAMP); 

•  Computational  meehanies  laboratories  at  UNIEEI; 

•  Computational  meehanies  laboratories  at  UNI  CAMP. 
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1.6  Air  Force  Relevance 

This  research  investigates  numerical  models  for  detection  of  damages  in  aircraft  structures.  The 
timely  detection  of  damage  is  an  important  feature  for  flight  safety.  This  research  is  also  relevant 
in  aircraft  maintenance,  with  the  possibility  for  using  such  on-board  structural  health  monitoring 
system  as  a  substitute  for  some  costly  structural  inspections  during  regular,  scheduled 
maintenance. 

1.7  Cost  /  Funding 

For  this  project,  a  grant  of  U$70,000.00  was  awarded  (U$  20,000.00  for  the  first  year,  and  U$ 
25,000.00  per  year,  for  the  second  and  third  years).  The  funding  was  directed  to  finance  research 
expenses  (both  for  the  researchers  and  the  students),  evenly  divided  for  each  university,  during 
these  three  years  of  collaborative  work.  Expenses  covered  included  stipends  for  students,  and 
also  general  expenses  related  to  the  project  (for  example,  travel  expenses,  computers,  software, 
computer  consumables,  books,  etc). 

1.8  Contact  Information  for  the  P.I. 

Prof.  Dr.  Ariosto  Bretanha  Jorge 
UNIFEI  -  Universidade  Eederal  de  Itajuba 
Institute  de  Engenharia  Mecanica 
Av.  BPS  1303 

37500-903  -  Itajuba,  MG  -  Brazil 

Tel  (Brazil):  +(35)  3629-1462  (office);  +(35)  3629-1152  (secretary); 

Pax:  +(35)  3629-1148  (secretary) 

Email:  ariosto.b.iorge@unifei.edu.br 

Home  page:  http://www.gemec.unifei.edu.br/ariosto/ 


2.  Research  description 

Modeling  the  damage  detection  problem:  general  aspects 

The  detection  of  damage  (or  failure)  in  structures  is  an  important  area  in  engineering,  with 
several  fields  of  application,  among  which  one  can  point  out:  flight  safety  and  aircraft 
maintenance,  piping  and  containers  in  the  oil  industry,  structures  in  nuclear  power  plant  industry, 
etc.  The  development  of  damage  detection  techniques  can  bring  up  technological  advances  in 
order  to  increase  the  structural  reliability  (safety),  contributing  to  a  better  structural  integrity 
analysis  and  to  a  better  evaluation  of  the  remaining  service  life  (or  useful  life)  of  a  structure.  The 
analysis  of  a  damaged  structure  must  involve  the  numerical  treatment  of  data  gathered  from 
sensors  spread  throughout  critical  points  in  the  structure,  and  the  comparison  of  this  data  with 
numerical  results  used  as  reference  (for  example,  results  from  the  same  structure  without  damage 
of  with  given  (known)  damage).  These  damage  detection  techniques  could  involve  monitoring 
(in  real  time)  of  the  integrity  of  structural  elements  which  are  critical  and  are  difficult  to  access 
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(poor  aeeessibility),  sueh  as  some  elements  in  an  airframe,  or  in  an  oil  piping,  or  in  nuelear  or  oil 
installations. 

2.1  Damage  detection:  overview  of  the  Direct  and  Inverse  Problems 

To  analyze  a  damage  deteetion  problem  in  a  strueture,  first  the  modeling  of  the  direet  problems 
is  required,  to  obtain  the  behavior  of  this  strueture  in  the  presenee  of  one  or  more  pre-established 
damages,  with  assumed  format  and  size,  and  at  given  positions  (see  referenees  [4]  e  [5]  for 
damage  deteetion  problems). 

In  this  work,  two  methods  of  analysis  were  given  partieular  attention,  for  the  direet  method: 

1)  The  study  of  the  stress  and  strain  distributions  in  struetural  elements,  performed  through  the 
BEM  (boundary  element  method)  modeling  of  elastostaties  problems  (for  defeets  modeled  as 
holes  in  the  strueture)  or  through  the  BEM  modeling  of  linear  elastie  fraeture  meehanies 
problems  (for  eraeked  struetures)  (see  referenees  [6]  to  [13]  for  the  boundary  element  method 
applied  to  fraeture  meehanies); 

2)  The  study  of  the  distribution  of  some  sealar  field  throughout  the  strueture,  modeled  using 
BEM.  Eor  problems  governed  by  the  Eaplaee  or  Poisson  equations,  a  potential  field,  sueh  as  the 
temperature  distribution,  is  the  quantity  of  interest.  Other  potential  fields,  sueh  as  the  sound 
pressure  in  the  strueture  due  to  sound  waves  (emitted  from  a  pre-established  souree)  were  also 
investigated.  The  presenee  of  the  defeet  in  the  strueture  infiuenees  the  distribution  of  these  sealar 
fields.  BEM  models  were  used  both  for  Eaplaee  /  Poisson  problems  (heat  transfer  by  eonduetion 
in  the  strueture)  and  for  the  modeling  of  the  aeoustie  problem  (see  referenees  [14]  to  [17]  for  the 
boundary  element  method  applied  to  aeousties). 

The  study  of  other  modeling  approaehes  is  being  planned  as  an  on-going  researeh  work,  and  is 
not  eovered  in  this  report.  The  methods  of  interest,  for  this  future  researeh,  inelude  damage 
deteetion  approaehes  based  on  vibration  analysis,  and  also  based  on  a  eombination  of  Einite 
Element  and  Boundary  Element  Methods  as  a  EEM/BEM  direet  model  (the  BEM  model  being 
best  suited  for  large  gradients  in  the  field  being  eonsidered,  while  the  EEM  aeeommodating  well 
properties/material  ehanges  throughout  the  strueture). 

Eor  the  study  of  the  inverse  problem,  the  model  eonsists  of  two  parts: 

a)  Monitoring  the  struetural  integrity,  from  experimental  measurements,  with  a  eertain 
number  of  sensors  spread  throughout  the  strueture.  With  this,  some  knowledge  is 
obtained  about  the  distribution  of  stresses  or  strains  (for  example,  by  means  of  strain- 
gages)  or  about  the  distribution  of  a  variable  derived  from  the  aeoustie  pressure 
throughout  the  strueture  (for  example,  by  means  of  mierophones,  aeeelerometers,  or  other 
sensors). 

b)  Computation  of  a  funetional  obtained  from  adding  differenees  (evaluated  in  all 
measurement  points)  between  the  values  evaluated  using  the  numerieal  model  from  the 
direet  problem,  for  a  given  damage  that  was  assumed,  and  the  experimental  values 
measured  in  the  same  points  for  the  strueture  with  the  real  damage  or  eraek.  This 
funetional  is  a  funetion  of  the  eraek  or  damage  loeation,  either  numerieal  or  measured 
from  the  real  strueture.  This  funetional  is  expeeted  to  inerease  in  value  when  the  assumed 
numerieal  defeet  is  far  away  from  the  real  defeet.  Also,  this  funetional  is  expeeted  to 
reaeh  its  minimum  value  when  both  defeets  (for  example,  the  numerieal  and  real  eraeks) 
eoineide.  Thus,  the  inverse  problem  is,  in  faet,  an  optimization  problem  for  the  seareh  of 
a  global  minimum  for  this  funetional.  In  this  work,  loeal  optimization  methods  (sueh  as 


7 


AFOSR/  SOARD.  Grant:  FA9550-06-1-0542.  Final  Report.  P.L:  Dr.  Ariosto  B  .Jorge. 

Numerieal  Modeling  of  Inverse  Problems  for  Damage  Deteetion  in  Aireraft  Struetures. 

linear  or  quadratie  sequential  programming)  were  used  to  eompare  with  global 
optimization  heuristies  (sueh  as  genetie  algorithms  and  differential  evolution). 

The  study  of  the  randomness  of  the  variables  involved  in  these  problems  is  of  great  importanee, 
in  order  to  the  eomputational  modeling  being  used  to  be  representative  of  the  real  strueture.  The 
study  of  the  randomness  was  ineipient  in  this  work,  using  stoehastie  optimization  teehniques.  To 
aeeount  for  these  uneertainties,  future  work  is  planned,  for  the  use  of  parameter  identifieation 
teehniques  (sueh  the  Kalman  filter  approaeh),  together  with  the  treatment  of  the  uneertainties  of 
the  variables  (using  fuzzy  optimization  teehniques)  or  the  treatment  of  the  randomness  of  these 
variables  (using  stoehastie  optimization  teehniques,  sueh  as  response  surfaees,  or  using  Monte 
Carlo  simulation). 

2.2  Direct  Problem 

The  boundary  element  method  is  a  numerieal  proeedure  well  adapted  for  the  modeling  of  a 
eraeked  strueture  made  of  an  isotropie  material  (aluminum,  for  example)  or  made  of  orthotropie 
materials  (sueh  as  eomposite  materials).  In  this  method,  the  distribution  of  the  quantities  of 
interest  in  the  domain  is  obtained  from  the  information  of  the  distribution  of  some  quantities  in 
the  boundary.  Thus,  in  this  method,  the  problem  is  deseribed  based  on  what  happens  in  its 
boundaries,  redueing  the  dimension  of  the  problem  and  simplifying  numerieally  the  treatment. 
Furthermore,  the  boundary  element  method  offers  an  additional  advantage  for  the  fraeture 
meehanies  problem,  as  the  sealar  and  veetor  fields  of  the  variables  of  interest  ean  be  deseribed 
with  reasonable  aeeuraey,  even  when  these  fields  are  singular,  as  it  is  the  ease  for  the  stress  field 
near  the  eraek  ends  (near  eraek  tip  in  2D,  or  near  eraek  eontour  in  3D). 

In  the  ease  of  the  approaeh  for  the  damage  deteetion  problem  made  by  means  of  the  analysis  of 
the  aeoustie  response  of  the  strueture  under  exeitation,  perturbations  in  the  expeeted  response 
imply  in  the  presenee  of  damage.  Thus,  the  defeets  and  damage  in  the  strueture  shall  eharaeterize 
its  dynamie  behavior.  The  modeling  of  the  direet  problem  may  inelude  also  the  study  of  the 
damage  evolution  through  time  (sueh  as,  for  example,  the  veloeity  of  eraek  propagation),  in 
order  to  estimate  the  remaining  useful  life  (safe  life)  of  the  strueture. 

2.3  Inverse  Problem 

The  inverse  problem  might  be  modeled  by  means  of  optimization  teehniques  or  by  identifieation 
teehniques.  In  the  following  diseussion,  some  aspeets  of  these  teehniques  are  detailed.  For  the 
diseussion,  a  simple  Laplaee  problem  for  the  distribution  of  a  potential  field  in  a  domain  is 
eonsidered.  The  damage  is  simulated  by  the  presenee  of  small  holes  in  the  domain,  and  the  goal 
is  to  obtain  the  size  (diameter  of  the  hole)  and  the  loeation  (veetor  position  of  the  eenter  of  the 
hole)  of  the  damage. 

The  direet  method  (BEM)  provides  one  pieee  of  information  (the  potential)  for  any  desired  point 
in  the  domain.  Without  the  hole,  the  distribution  of  the  potential  is  known  a  priori.  If  a  small  hole 
is  ineluded,  the  potential  distribution  is  unknown  and  must  be  obtained  numerieally  from  the 
BEM  solution.  The  goal  in  this  problem  is  to  implement  two  inverse  methods  (optimization  and 
identifieation)  and  to  diseuss  the  diffieulties  in  the  implementation  and  advantages  of  eaeh 
method,  to  find  out  whieh  one  is  more  appropriate  to  solve  the  problem.  The  results  obtained  for 
the  inverse  method  by  means  of  this  two  independent  teehniques  (genetie  algorithm  (GA)  for  the 
optimization  proeedure  and  artifieial  neural  networks  (ANN’s)  for  the  identifieation  proeedure) 
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are  eompared  to  analyze  the  effieieney  of  eaeh  method  in  finding  the  loeation  and  dimension  of 
the  hole. 

Inereasing  the  problem  eomplexity,  the  BEM  for  the  elastieity  problem  ean  be  used,  with  some 
given  boundary  eonditions  for  the  displaeement  and  traetion.  Differently  from  the  BEM  for  the 
potential,  the  BEM  for  elastieity  (in  a  2D  problem)  provides  two  pieees  of  information  at  a  single 
interior  point  -  one  normal  stress  and  one  shear  stress.  But  this  information  eannot  be  used 
direetly  in  the  optimization  problem,  as  it  depends  on  the  system  of  eoordinates  being  used,  or  on 
the  normal  direetion  of  the  eutting  plane  that  passes  through  the  point  of  interest.  Therefore,  a 
ehoiee  was  made  to  adopt  the  stress  invariants  of  the  stress  tensor  at  the  point  of  interest-  in  2D, 
the  mean  stress  and  the  oetahedral  stress  -  as  the  veetor  field  to  be  analyzed  and  used  in  the 
optimization  problem. 

A  eomparison  of  different  (and  independent)  optimization  and  identifieation  teehniques  for  this 
inverse  problem  of  damage  deteetion  is  also  desirable,  to  eheek  for  the  robustness  of  the  different 
approaehes  in  finding  the  damage.  The  global  optimization  teehniques  that  ean  be  used  inelude: 
genetie  algorithm  (GA),  differential  evolution  (DE),  ant  eolony  (AC),  ete.  The  parameter 
identifieation  teehniques  that  ean  be  used  inelude:  artifieial  neural  networks  (ANN’s)  and 
Kalman  filter  (KE). 

In  order  to  inerease  the  eomplexity  of  the  damage  deteetion  problem  using  the  BEM  formulation 
for  elastieity,  the  goal  eould  be  to  identify  and  loeate  (using  one  or  more  of  the  inverse  problems 
approaehes)  the  presenee  in  the  plate  of  one  or  more  eireular  holes  (number  of  holes,  radius  and 
loeation  of  eaeh  hole)  and  also  one  or  more  ellipses  (number  of  ellipses,  axis  orientation  and  size 
(large  and  small  axis)  and  loeation  of  eenter  for  eaeh  ellipse).  To  inerease  even  further  the 
eomplexity  of  the  damage  deteetion  modeling,  the  BEM  for  fraeture  meehanies  ean  be  used,  in 
order  to  loeate  eraeks  in  the  plate  (number  of  eraeks,  and  their  size,  orientation,  and  loeation). 

Eor  all  the  different  direet  problems,  the  same  global  optimization  teehniques  and/or  parameter 
identifieation  proeedures  ean  be  used  to  solve  the  inverse  problem  for  damage  deteetion.  Erom 
the  point  of  view  of  the  inverse  problem,  the  direet  model  is  just  a  ‘black  box'  to  be  supplied  to 
give  the  numerieal  information  needed  to  be  used  in  the  optimization  or  identifieation  proeedure. 
The  starting  point  in  the  researeh  is  to  analyze  and  diseuss  the  deteetion  of  just  one  hole  in  the 
strueture.  Eater  in  the  researeh,  by  modifying  the  implementation  of  the  inverse  problem,  the 
model  would  be  able  to  deteet  more  than  one  damage  in  a  partieular  struetural  element. 
Coneomitant  to  this  researeh,  several  different  eonfigurations  for  the  direet  problem  are 
evaluated  (BEM  for  eraeked  orthotropie  materials,  material  interfaees,  ete).  The  goal  is  to 
implement  BEM  eodes  representing  situations  that  emulate  better  a  real  aireraft  strueture 
(ineluding  struetural  patehes,  reinforeements,  ete).  These  BEM  formulations  for  the  direet 
problem  are  being  developed  mostly  by  the  researehers  at  UNICAMP.  On  the  other  hand,  the 
implementation  of  the  inverse  method  is  being  developed  mostly  at  UNIEEI.  The  two  researeh 
fronts  are  under  way  in  parallel,  so  that  the  final  eode  would  inelude  as  many  direet  eodes  as 
possible  (to  allow  for  a  representative  number  of  possible  eombinations  of  damages  /  airframes 
to  be  evaluated),  as  well  as  a  number  of  different  (independent)  inverse  models,  to  inerease  the 
reliability  in  the  damage  information  (quantity,  loeation,  and  size). 

Damage  detection  by  means  of  optimization  techniques 

In  an  experimental  analysis,  the  data  gathered  eome  from  sensors  spread  throughout  the 
strueture,  loeated  at  a  number  of  points.  The  experimental  analysis  is  not  being  undertaken  at  this 
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moment,  thus  the  size  and  location  of  the  real  damage  in  the  plate  are  being  assumed  and 
simulated  by  using  the  BEM.  In  this  first  example  problem,  BEM  for  the  potential  is  being  used, 
and  the  temperature  values  at  some  interior  points  replace  the  information  that  would  have  been 
collected  by  the  sensors  at  these  points  in  the  plate. 

In  order  to  solve  the  inverse  problem  for  damage  detection,  an  optimization  algorithm  (GA,  for 
example)  is  used.  The  evaluation  function  (fitness  function)  for  the  GA  is  formulated  as  a 
functional  defined  as  a  difference  between  measured  values  (simulated,  in  this  case)  of  the  local 
difference  in  the  potential  (between  the  undamaged  plate  and  the  plate  with  the  damage)  and  the 
values  of  the  same  differences  in  potential  calculated  at  the  same  points  by  the  damage  detection 
code  (assuming  several  different  locations  and  sizes  for  the  ‘numerical’  damage).  The 
minimization  of  the  functional  allows  the  damage  detection  code  to  find  the  unknown  parameters 
of  the  damage.  A  general  formulation  for  the  functional  is  shown  in  Equation  (1). 


1  ” 

J j  =  —  '^(measured I  -  calculated 

2  ,=i 


(1) 


Where: 

n  -  Number  of  sensors  placed  in  the  plate  (number  of  internal  points  i  where  differences 
are  evaluated); 

measuredi  -  Vector  of  simulated  values  for  the  differences  in  potential  (obtained  using 
BEM,  these  values  represent  the  values  measured  in  the  plate  points  for  a  given  damage  location 
and  size); 

calculatedji  -  Vector  of  values  for  the  differences  in  potential  calculated  by  the  damage 
detection  code  for  each  individual  j. 

Eigure  1  represents  an  undamaged  thin  plate  with  four  sensors  indicating  the  points  where  the 
measurement  of  the  quantities  of  interest  (differences  in  potential,  in  this  case,  or  stresses,  for 
example,  in  the  elasticity  case)  is  being  performed. 


Figure  1  -  Undamaged  plate  with  four  representative  sensors 

In  order  to  solve  the  damage  detection  problem,  an  initial  population  is  given  to  the  GA.  This 
initial  population  is  formed  by  individuals  which  constitute  a  possible  solution  for  the  problem. 
These  individuals  are  represented  by  chromosomes  which  are  themselves  constituted  by  genes. 
Each  gene  in  a  chromosome  represents  one  variable  in  the  problem  (for  example,  the  x  and  y 
coordinates  and  the  radius  r  of  the  hole).  As  an  example,  Eigure  2(a)  to  2(c)  represents  three 
possible  configurations  of  chromosomes.  While  the  location  and  size  of  the  hole  varies,  the 
number  of  sensors  and  their  locations  are  always  the  same,  for  all  chromosomes  (these  are  also 
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the  locations  of  the  ‘real’  sensors  for  the  measurements).  The  information  on  the  quantity  of 
interest  is  collected  at  these  sensor  locations  for  all  cases. 


(a)  chromosome  1  (b)  chromosome  2  (c)  chromosome  3 

Figure  2  -  Plate  with  a  hole:  three  possible  eonfigurations  for  the  ehromosomes 


Some  of  the  plate  properties  (such  as  material  properties  and  geometry)  may  not  be  known 
exactly  a  priori  (due  to  randomness  in  the  manufacturing  or  fabrication  processes,  for  example). 
Material  properties  (such  as  elasticity  modulus  E,  and  Poisson’s  ratio  v,  in  the  elastic  problem) 
and  geometric  parameters  (see  Figure  3,  for  example)  are  suited  to  be  obtained  in  the  real 
structure  by  means  of  parameter  identification  procedures  (such  as  ANN,  for  example). 


Figure  3  -  Plate  geometry:  possible  parameters  to  be  identified 

Other  variables  (such  as  the  loading  (shown  in  Figure  4)  and  the  boundary  conditions)  may 
contain  uncertainties  and  randomness  and  may  need  to  be  treated  as  random  variables  in  the 
damage  detection  code. 


loading 


Figure  4  -  Loading:  possible  variables  to  be  treated  as  random 
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Also,  the  number,  size  and  location  of  the  different  damages  may  need  to  be  treated  either  by 
parameter  identification  or  as  random  variables. 

For  a  straight  crack,  its  parameters  may  be  (see  Figure  5): 

-  size  =  a; 

-  orientation  =  6] 

-  position  =  (xo,yo)- 


Figure  5  -  Crack  parameters:  size,  orientation,  position 

The  treatment  of  randomness  is  detailed  in  the  next  session  (Modeling  of  Uncertainties).  A  result 
expected  from  the  on-going  research  is  the  discussion  on  which  method  for  modeling  of 
uncertainties  (identification  procedures  versus  probabilistic  methods)  is  best  suited  for  each  of 
the  different  variables  and  parameters  of  the  problems. 

For  the  transition  from  the  plate  with  a  hole  to  a  cracked  plate,  the  direct  model  using  BEM  is  the 
only  subroutine  that  changes  in  the  code  (from  BEM  for  elasticity  to  BEM  for  linear  elastic 
fracture  mechanics).  The  procedures  for  the  inverse  problem  (either  optimization  or 
identification  techniques)  remain  unchanged.  Thus,  for  the  simulation,  the  position,  orientation 
and  size  of  the  ‘real’  crack  are  assumed  (see  Figure  6).  Also,  the  numerical  solutions  for  the 
different  chromosomes  are  obtained  in  the  same  way  as  before  (see  Figure  7  for  three  possible 
chromosomes). 


Figure  6  -  Cracked  plate:  size,  orientation,  position  of  the  ‘real’  crack  (simulated,  in  this  case) 
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2.4  Modeling  of  Uncertainties 

This  part  of  the  research  involves  the  numerical  modeling  of  an  engineering  optimization 
problem  with  two  basic  features: 

i)  the  various  parameters  and  variables  of  the  system  being  studied  are  not  deterministic, 
and  a  proper  treatment  of  this  variability  leads  to  robust  optimization  techniques,  where 
the  goal  is  to  obtain  not  only  optimum  values  for  the  objective  functions,  but  also 
minimum  variations  in  these  objective  functions  in  the  neighborhood  of  these  optimum 
points.  In  this  context,  the  words  'stochastic  optimization'  and  'robust  optimization'  lead 
to  the  same  idea,  as  the  goal  of  the  treatment  of  the  stochasticity  of  the  variables  and 
parameters  of  the  problem  is  to  obtain  robust  optima.  In  this  case,  the  optima  are  points  in 
the  feasible  region,  wherein  the  values  of  the  objective  functions  are  insensitive  to 
variations  around  these  points; 

ii)  the  fact  that  the  model  has  to  look  not  only  for  the  optimum  values  of  the  objective 
functions  but  also  for  robustness  (that  is,  a  small  variability  of  these  objective  functions 
around  these  optima)  shows  that,  in  each  problem,  there  is  always  more  than  one 
objective  function.  Thus,  there  is  a  need  for  decision-making  procedures  with  respect  to 
these  multiple  objectives,  which  may  involve  the  use  of  different  multiple-objective 
optimization  techniques,  such  as  weighting,  prioritizing  (goal  programming),  the  use  of 
objective  functions  as  constraint  equations,  the  use  of  fuzzy  membership  functions, 
obtaining  Pareto  limiting  regions  or  curves,  etc. 

Thus,  the  modeling  of  uncertainties  of  the  damage  detection  problem  in  an  aeronautical  structure 
involves  stochastic  multi-objective  optimization  techniques  in  the  modeling  of  the  inverse 
problem. 

Modeling  the  multiple-objective  optimization  problem:  general  aspects 

The  traditional  optimization  methods  usually  treat  the  variables  of  the  problems  as  deterministic. 
For  a  review  of  traditional  calculus-based  algorithms  (for  the  search  of  local  optima),  see 
references  [18]  and  [19]. 

Heuristics  that  search  for  global  optima  have  been  proposed  in  the  literature,  several  of  which 
based  in  the  imitation  of  behaviors  found  in  nature.  An  example  of  this  is  the  'survival  of  the 
fittest',  found  in  heuristics  such  as  evolutionary  algorithms,  genetic  algorithms,  differential 
evolution,  particle  swarm  optimization,  etc  (see  references  [20]  to  [27]). 
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But  in  several  eases,  these  algorithms  still  eonsider  as  deterministie  the  problem  variables  (for 
example,  the  loading,  boundary  eonditions,  material  properties,  geometry  ete)  and/or  parameters 
(for  example,  the  eoefficients  in  the  objeetive  funetions  or  in  the  eonstraint  equations). 

Treatment  of  stochasticity  of  variables  and  parameters 

The  modeling  of  stoehastieity  both  in  the  variables  and  in  the  parameters  of  a  problem  involves 
the  use  of  probabilistie  methods  in  engineering,  sueh  as  Monte  Carlo  simulation.  Response 
Surfaee  teehniques.  Design  of  Experiments,  First  Order  Reliability  Methods  (FORM)  or  Seeond 
Order  Reliability  Methods  (SORM),  logistie  regression,  ete. 

The  stoehastieity  ean  be  also  used  in  identifieation  proeedures  in  two  steps:  first,  a  set  of  random 
information  is  used  to  identify  the  system  parameters  (for  example,  via  Kalman  filter),  and 
seeond,  an  independent  set  of  stoehastie  data  is  used  for  the  treatment  of  the  random  variables  of 
the  problem  (see  referenees  [28]  to  [32]  on  different  probabilistie  methods  in  engineering). 

Decision  techniques  in  the  treatment  of  multiple  objectives 

Decision  techniques  regarding  the  multiple  objectives  of  an  optimization  problem  may  include: 
weighting,  assignment  of  priorities,  the  use  of  objective  functions  as  constraint  equations,  the  use 
of  fuzzy  membership  functions  in  the  decision-making  process,  obtaining  regions  or  curves  of 
Pareto  limits,  etc. 

In  certain  cases,  objectives  of  different  natures  may  need  to  be  considered,  and  their  combination 
(through  weighting,  for  example)  may  not  be  possible.  For  example,  on  a  particular  problem,  one 
objective  may  happen  to  be  written  as  a  real  function,  while  another  objective  may  involve  only 
integer  numbers,  and  a  third  objective  may  involve  only  a  qualitative  response.  In  this  case,  a 
promising  technique  for  a  proper  combination  of  these  objectives  of  different  natures  could  be 
the  use  of  fuzzy  membership  functions  for  each  objective  function,  looking  for  the  optimization 
of  one  function  only,  namely,  the  summation  of  all  the  fuzzy  membership  functions.  This 
technique  could  even  allow  the  designer  to  include  a  bias  through  one  or  another  objective 
function,  if  this  is  considered  necessary  (see  references  [33]  and  [34]  on  the  use  of  fuzzy  logic  in 
optimization). 

3.  Main  accomplishments 

This  research  is  a  collaborative  effort  between  the  Computational  Mechanics  groups  at  UNIFEI 
(Itajuba)  and  UNICAMP  (Campinas).  The  research  on  the  direct  problem  has  been  performed 
under  supervision  of  Professors  Paulo  Sollero  and  Eder  Eima,  at  UNI  CAMP,  while  the  research 
on  the  inverse  problem  and  on  the  modeling  of  uncertainties  has  been  performed  under 
supervision  of  Professors  Ariosto  and  Sebastiao  Simoes,  atUNIEEI. 

As  results  directly  related  to  this  research,  several  publications  and  monographs  were  obtained, 
as  detailed  in  the  Appendix.  The  publications  were  concentrated  on  conference  papers  and 
journal  articles,  while  the  monographs  were  concentrated  on  thesis  and  dissertations  defended  by 
the  students  working  on  the  research  groups,  both  at  UNIEEI  and  UNICAMP.  Besides  the 
published  work  and  the  defended  thesis  and  dissertations,  some  research  papers  and  student 
thesis/dissertation  defenses  are  also  expected  to  occur  in  the  near  future,  related  to  this  work.  The 
list  in  the  Appendix  includes  the  on-going  research,  leading  to  student  dissertations  /  thesis. 
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3.1  Main  accomplishments  on  the  direct  and  inverse  models 

Throughout  the  research  project,  the  complexity  of  the  modeling  of  the  direct  and  inverse 
problems  has  increased,  pursuing  two  goals: 

•  The  goal  for  the  direct  model  is  to  be  improved  in  order  to  reproduce  as  close  as  possible  the 
reality  of  an  aircraft  structure. 

•  The  goal  for  the  inverse  model  is  to  be  as  reliable  as  possible.  For  that,  a  combination  of 
independent  optimization  and  identification  procedures  is  deemed  necessary.  The  reliability 
of  the  damage  detection  procedure  is  expected  to  be  high,  if  two  or  more  independent 
approaches  indicate  the  same  location  and  size  for  the  structural  damage. 

In  addition,  the  current  inverse  models  already  include  some  stochastic  modeling,  but  the  proper 
treatment  of  the  uncertainties  is  an  on-going  research,  and  is  also  the  object  of  collaborative 
research  planned  for  the  near  future. 

The  main  accomplishments  of  this  research  project  can  be  summarized  as  follows: 

•  For  the  direct  problem: 

o  BEM  models  implemented:  cracked  composite  plates  with  repair  patches  (static  and 
dynamic),  evaluation  of  adhesive  shear  stress,  anisotropic  fundamental  solution  using 
numerical  integration  (Radon  transformation); 
o  Comparison  between  Dual  Reciprocity  BEM  and  cell  domain  integration  to  treat 
remaining  domain  integrals  (due  to  the  shear  interaction  forces). 

•  Eor  the  inverse  problem: 

o  Damage  detection  algorithm  implemented  for  standard  BEM  models  (potential, 
elasticity,  acoustics); 

o  Optimization  methods  (Sequential  Quadratic  Programming  -  SQP;  Genetic 
Algorithms  -  GA)  and  identification  techniques  (Artificial  Neural  Networks  -  ANN) 
were  compared  for  structures  with  deterministic  parameters. 

The  models  implemented  are  detailed  in  the  several  published  works  (conference  papers  and 
journal  articles)  cited  in  the  Appendix.  Some  illustrative  results  for  the  direct  model  research  are 
shown  below,  in  Eigures  8  to  11.  Eigure  8  shows  the  boundary  element  model  for  a  cracked  plate 
with  a  composite  patch.  The  remaining  domain  integrals  (due  to  the  shear  interaction  forces) 
need  to  be  evaluated  either  by  a  Dual  Reciprocity  approach  of  by  cell  domain  integration. 


Figure  8:  Circular  composite  patch  over  a  cracked  square  sheet:  BEM  model 
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Figure  9  shows  the  stress  distribution  obtained  with  the  BEM  model.  The  presence  of  the  circular 
patch  has  alleviated  the  stresses  on  the  cracked  region  of  the  plate. 


Figure  9:  Circular  composite  patch  over  a  cracked  square  sheet:  stress  results  using  BEM 


Eigure  10  shows  a  representation  of  the  problem  being  modeled,  in  which  there  is  an  adhesive 
layer  between  the  composite  patch  (anisotropic)  and  the  metallic  plate  (isotropic).  The  adhesive 
layer  was  treated  in  the  BEM  model  implemented. 


BB  Isotropic  plate 

■  Anisotropoic  repair 

■  Adhesive  layer 

Figure  10:  Dynamic  analysis  model:  plate,  composite  patch,  adhesive  layer 


Section  A-A 


Figure  1 1  shows  the  numerical  integration  results  to  obtain  the  anisotropic  fundamental  solution 
in  3D  using  the  Radon  transformation,  pointing  out  the  need  to  increase  the  number  of  Gauss 
points  in  this  numerical  integration,  in  order  to  obtain  a  proper  reconstruction  of  the  smooth 
anisotropic  fundamental  solution,  to  be  used  with  the  BEM  formulation. 
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4  Gauss  points  30  Gauss  points 


Figure  11:  Anisotropic  fundamental  solution  3D:  numerical  integration  requires  30  or  more 
Gauss  points. 

Regarding  the  research  on  the  inverse  model,  some  illustrative  results  are  shown  below,  in 
Figures  12  to  14.  Figure  12  shows  the  convergence  pattern  of  the  inverse  method  using  the 
genetic  algorithm  (GA)  as  the  optimization  approach,  and  the  BEM  model  of  the  Laplace 
equation  as  the  direct  model  for  the  potential  distribution  along  the  plate.  Besides  the  usual 
operations  of  mutation  and  cross-over,  the  use  of  elitism  improves  the  accuracy  of  the  damage 
localization  results. 


Figure  12:  Influence  of  GA  parameters:  high  elitism  is  better 

Figure  13  shows  the  convergence  pattern  of  the  inverse  method  using  the  genetic  algorithm  (GA) 
as  the  optimization  approach,  and  the  elastostatics  BEM  as  the  direct  model.  The  localization 
results  present  smaller  variability  when  the  distribution  of  a  scalar  quantity  along  the  plate  is 
used  (in  this  case,  the  mean  stress  or  the  octahedral  stress,  which  are  the  invariants  of  the  stress 
tensor),  instead  of  a  vector  quantity  (in  this  case,  the  stress  components,  which  are  direction- 
dependent). 
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Figure  13:  Stress  invariants:  better  results  (less  variability)  than  stress  veetor  eomponents 

For  the  aeousties  model,  the  aeoustie  pressure  is  a  potential  field.  This  aeoustie  potential  at  a 
partieular  loeation  (the  sensor  loeation)  is  a  funetion  of  time,  and  ean  be  measured  for  the  ease  of 
the  real  hole,  and  ean  be  simulated  using  BEM,  for  the  ease  of  the  numerieal  hole.  A  funetional 
ean  be  built  as  differenees  in  the  areas  below  the  eurves  of  the  potential  for  the  real  and  measured 
eases.  Figure  14  shows  that  this  funetional  eorrelates  direetly  with  the  distanees  between  the 
numerieal  and  the  real  hole  (the  value  of  this  funetional  gets  smaller  as  the  numerieal  hole 
approaehes  the  real  hole). 


Figure  14:  Aeousties  model:  “numerieal”  holes  approaehing  the  “reaf  ’  hole.  The  inverse  model 
identifies  the  numerieal  hole  elosest  to  real  hole. 

3.2  Establishment  of  collaborative  research  work  (on-going  work) 

The  researeh  on  the  modeling  and  simulation  of  the  inverse  problem  has  involved  some 
eollaborative  work  with  Prof.  G.  Walker  and  with  Professor  P.  K.  Basu,  both  from  Vanderbilt 
University  (Nashville,  TN).  This  eollaboration  was  very  important  for  diseussing  modeling 
teehniques  for  the  inverse  problem.  Furthermore,  this  eollaboration  has  ereated  a  synergy 
between  the  eomputational  meehanies  group  from  UNIFEI  in  Brazil  and  Vanderbilt  University, 
important  for  the  planned  modeling  of  the  uneertainties  of  the  damage  deteetion  problem. 
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To  aoeomplish  this  eollaborative  work,  Prof  Ariosto  spent  some  time  in  2009  at  Vanderbilt 
University,  as  a  visiting  scholar,  while  continuing  to  coordinate  the  research  work  being  done  by 
the  students  and  researchers  from  the  UNIFEI  research  group,  in  Brazil.  Funding  for  the 
expenses  of  this  visit  to  Vanderbilt  University  was  obtained  through  a  Brazilian  agency,  CNPq. 
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Abstract.  A  boundary  element  formulation  for  the  analysis  of  isotropic  cracked  sheets,  repaired 
with  adhesively  bonded  anisotropic  patches  is  presented.  The  sheet  and  the  patch  are  modeled  using 
the  boundary  element  method.  The  crack  in  the  isotropic  sheet  is  modeled  using  the  dual  boundary 
element  method.  The  interaction  between  the  isotropic  sheet  and  the  patch  is  modeled  considering 
shear  body  forces  uniformly  distributed  on  the  interaction  zone  using  a  linear  elastic  relationship.  Two 
different  techniques  are  used  in  the  present  boundary  element  implementation  to  treat  the  domain 
integrals  that  arise  in  the  formulation  due  to  shear  interaction  forces.  These  techniques  are  the  cell 
domain  integration  and  the  dual  reciprocity  boundary  element  method.  Examples  show  that  results 
obtained  for  the  shear  stress  distribution  in  the  adhesive  layer  are  in  good  agreement  with  analytical 
solution. 

Introduction 

Adhesively  bonded  composite  patches  are  increasingly  used  in  aircraft  structure  repairs  in  order  to 
extend  the  life  of  cracked  structures  and  avoid  high  expenses  owing  to  the  replacement  of  cracked 
components.  In  aeronautical  applications,  when  a  non-destructive  technique  detects  a  crack,  it  is 
usually  necessary  to  drill  the  crack  tip  region  in  order  to  decrease  the  stress  concentration  and  then 
apply  a  layer  of  adhesive  patch  on  this  region  to  avoid  the  crack  growth.  The  patch  transfers  the  load 
from  the  cracked  structure  to  the  repair,  avoiding  crack  opening  and  crack  propagation.  The  main 
advantages  of  bonded  patches,  when  compared  to  other  types  of  repairs  such  as  riveted  patches,  are 
the  homogeneous  load  transfer  between  the  cracked  plate  and  the  repair  and  the  absence  of  holes, 
which  are  stress  concentrators,  as  shown  by  Rose  and  Wang  [1]. 

Bonded  patches  in  cracked  structures  have  been  studied  by  many  researchers.  In  general,  the  sheet, 
the  patch,  and  the  adhesive  layer  are  considered  to  be  thin,  so  that  the  whole  component  does  not 
bend  out  of  its  plane,  and  the  problem  can  be  solved  using  the  two  dimensional  elasticity  theory.  The 
initial  works  analyzing  isotropic  patches  in  structures  were  presented  by  Erdogan  and  Arin  [2]  and 
Ratwani  [3],  in  the  seventies.  These  works  presented  the  study  of  bonded  repairs  in  infinite  plates  with 
cracks.  They  used  analytical  solutions  for  the  deformations  and  displacement  compatibility  between 
the  cracked  plate  and  the  repair. 

Mitchell,  Wooley  and  Chwiruth  [4]  used  the  finite  element  method  (EEM)  to  study  the  reinforce¬ 
ment  of  plates  induced  by  the  application  of  repairs.  They  used  two-dimensional  finite  elements  with 
constant  stress  distribution  and  the  plate  and  repair  were  coupled  through  nodes  where  conditions  of 
displacement  compatibility  were  imposed.  They  also  analyzed  the  presence  of  a  crack  in  the  plate. 
However,  they  did  not  consider  the  stress  singularity  at  the  crack  tip  and  did  not  evaluate  the  stress 
intensity  factors.  Jones  and  Callinan  (see  References  [5],  [6],  [7])  used  the  EEM  for  the  analysis  of 
metallic  plates  repaired  with  a  layer  of  composite  material.  They  developed  a  stiffness  matrix  to 
couple  the  plate,  the  adhesive  layer,  and  the  composite  repair.  Special  singular  elements  were  used  at 
the  crack  tip. 

Young,  Cartwright  and  Rooke  [8]  modeled  the  cracked  plate  and  the  repair  using  the  boundary 
element  method  (BEM).  Shear  stresses  in  the  adhesive  layer  and  body  forces  acting  on  the  plate  and 
on  the  repair  were  modeled  using  the  cell  integration  technique.  A  special  Green  function  for  domains 
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Figure  1:  Cracked  sheet  repaired  with  adhesive  patch 


with  cracks  was  used  to  model  straight  cracks,  which  limits  the  applicability  of  the  model. 

Tarn  and  Shek  [9]  studied  the  problem  of  cracked  plates  repaired  with  bonded  composite  patches. 
A  spring  was  nsed  to  conple  the  cracked  plate  model  with  the  repair  model.  The  repair  was  modeled 
using  the  FEM  and  the  crack  using  the  BEM.  Young  [10]  modeled  the  distributed  interaction  force 
between  the  plate  and  the  repair  by  discretizing  the  bonded  repaired  area  nsing  internal  cells  in  the 
boundary  element  formulation. 

Salgado  and  Aliabadi  [12]  introdnced  the  dnal  bonndary  element  method  (DBEM)  to  model  the 
metallic  cracked  plate  and  the  boundary  element  method  to  model  the  metallic  repair.  The  distributed 
forces  between  the  plate  and  the  repair  were  modeled  using  the  dual  reciprocity  boundary  element 
method  (DRBEM).  This  formnlation  was  applied  by  Salgado  and  Aliabadi  [13]  to  the  analysis  of 
metallic  thin  plates  reinforced  with  bonded  isotropic  repairs.  The  reinforced  plate  was  modeled  using 
the  BEM.  Shear  stresses  in  the  adhesive  layer  were  modeled  as  action-reaction  body  forces  exchanged 
between  the  plate  and  the  repair.  Widagdo  and  Aliabadi  [14]  extended  this  formulation  to  model  me¬ 
chanically  fastened  composite  repair  patches.  The  fasteners  were  considered  as  linear  springs  conpling 
the  cracked  sheet  and  the  anisotropic  repair  and  their  interaction  loading  was  modeled  as  a  summation 
of  discrete  point  forces.  Widagdo  and  Aliabadi  [15]  apply  this  formnlation  for  the  analysis  of  cracked 
sheet  repaired  with  adhesively  bonded  orthotropic  repairs. 

The  cnrrent  work  analyses  a  composite  repair  patch  adhesively  bonded  in  a  metallic  cracked  sheet. 
The  DBEM  is  used  to  model  the  isotropic  cracked  sheet  and  the  BEM  is  used  to  model  the  anisotropic 
composite  patch.  The  interaction  loading  between  the  sheet  and  the  patch  is  modeled  considering 
the  shear  forces  in  the  adhesive  layer  nniformly  distribnted  nsing  a  linear  elastic  relationship.  Two 
different  techniques  are  used  to  treat  domain  the  integrals  that  arise  in  the  formulation  due  to  the 
interaction  shear  forces:  the  cell  domain  integration  and  the  DRBEM. 

Numerical  examples  of  the  adhesive  stress  analysis  in  cracked  plate,  repaired  with  a  circular  and 
rectangnlar  composite  patches,  are  presented.  The  shear  stress  distribntions  obtained  with  the  cnrrent 
techniques  are  compared  to  the  analytical  solution  of  Rose  [16]  with  good  agreement.  Stress  intensity 
factors  are  calcnled  nsing  the  displacement  extrapolation  techniqne. 

1  Boundary  element  formulation 

Figure  1  presents  a  finite  isotropic  sheet,  containing  an  inner  crack  and  an  adhesive  patch.  In  this  case, 
the  interaction  forces  can  be  treated  as  nnknown  body  forces  exchanged  by  the  sheet  and  the  patch  in 
the  attachment  sub-region.  Considering  that  the  sheet  and  the  patch  remain  fiat  after  deformation, 
the  two-dimensional  elasticity  theory  can  be  nsed  to  model  this  problem.  In  this  case,  displacements 
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at  the  sheet  and  at  the  patch  have  to  be  compatible  with  the  shear  deformation  of  the  adhesive  layer 
connecting  them. 

When  the  sheet  is  deformed  due  to  applied  loads  on  its  boundaries,  interaction  forces  occur  be¬ 
tween  the  sheet,  with  contour  and  the  repair  patch,  with  contour  (see  figure  1).  In  this 
two-dimensional  case,  interaction  forces  in  the  plate  directly  underneath  the  repair  patch,  and  in  the 
patch  itself,  can  be  treated  as  unknown  body  forces  (action-reaction  pair).  As  shown  by  Salgado  and 
Aliabadi  [13],  the  boundary  integral  equation  for  the  displacement  of  a  source  point  x’  on  the  sheet  is 
given  by: 

cfj  (x')  Uj  (x')  -I-  j  T*f  (x',  x)  Uj  (x)  dT  =  j  U*f  (x',  x)  tj  (x)  dV  + 

Ts  Ts 

J  u*f  {x',x)bj  {x)dnE  i,j  =  1,2  (1) 


where  cfj  is  a  coefficient  which  depends  on  the  position  of  the  source  point  in  relation  to  the  boundary  of 
the  sheet  F5;  {x.' ,-x.)  and  Tfj^{x',x)  are  Kelvin’s  isotropic  fundamental  solutions  for  displacements 

and  tractions,  respectively,  for  the  two-dimensional  sheet  media;  uf  and  tf  are  displacement  and 
traction  vectors  at  the  boundary  of  the  sheet;  bf  are  interaction  forces  exchanged  between  the  sheet 
and  the  patch  in  the  domain  Qr  of  the  patch;  hs  is  the  thickness  of  the  sheet. 

Similarly,  the  displacement  of  a  source  point  x'  on  the  repair  is  given  by: 

cfj  (x')  uf  (x')  -I-  j  T*j^  (x',  x)  uf  (x')  dV  =  j  U*^  (x',  x)  tf  (x')  dF  -|- 
r  r 

^  2  (2) 


where  cfj  is  a  coefficient  which  depends  on  the  position  of  the  source  point  in  relation  to  the  boundary  of 
the  sheet  F^;  Uf^{x',x)  and  Tj*^(x',x)  are  anisotropic  fundamental  solutions  for  the  two-dimensional 
composite  repair;  uf  and  tf  are  displacement  and  traction  vectors  at  the  boundary  of  the  repair;  bf 
are  the  interaction  forces  exchanged  between  the  sheet  and  the  patch  in  the  domain  Qr  of  the  patch; 
Hr  is  the  thickness  of  the  sheet. 

In  this  work,  the  anisotropic  fundamental  solutions  for  two-dimensional  elastic  media  was  used  to 
model  the  mechanical  response  of  the  composite  patch  (see  Aliabadi  and  Sollero  [17]). 

The  crack  in  the  isotropic  sheet  was  modeled  using  the  DEEM.  The  traction  integral  equation  is 
applied  in  one  of  the  crack  faces  and  the  displacement  integral  equation  is  applied  in  the  other  crack 
face.  The  traction  integral  equation  is  given  by: 

i^')  +  Ri  i^')  J  Sfjk  (x',  x)uf  (x)  dr  =  Hi  (x')  J  D*fp.  (x',  x)  tf  (x)  dF 

Ts  Ts 

^  J  D*fk  (x',  x)  bf  (x)  dnR  i,j  =  1,2  (3) 

^R 

where  S*fj.{x',x)  and  are  linear  combinations  of  derivatives  of  fundamentals  solutions  for 

traction  and  displacement  Tj*^(x',  x)  and  x),  respectively,  and  n*  are  the  components  of  a  unit 

vector  outward  to  the  boundary  in  the  collocation  point. 

Now,  considering  a  uniform  shear  deformation  through  the  adhesive  thickness,  as  proposed  by 
Salgado  and  Aliabadi  [13],  and  neglecting  shear  deformations  in  the  sheet  and  in  the  patch,  the  body 
force  bj{x'),  that  is  equal  to  the  shear  stress  in  the  adhesive  Tj(x'),  can  be  written  as  a  function  of 
the  difference  Auj  between  the  displacements  uf  of  a  point  x'  (x'  G  Qr)  on  the  sheet  and  uf  of  a 
corresponding  point  on  the  repair  patch,  as: 
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(x')  =  Tj  (x') 


J  =  l,2 


(4) 


where  is  the  thickness  of  the  adhesive  layer,  Ga  is  the  transversal  stiffness  modnlns  of  the  adhesive 
material. 


2  Domain  integral  techniques 

As  can  be  seen,  eqnations  (1)  and  (2)  reqnire  the  calcnlation  of  domain  integrals.  Two  different 
techniques  were  used  and  compared  to  treat  the  domain  integrals  that  arise  in  the  formulation  due 
to  the  shear  interaction  forces.  These  techniques  are  the  cell  domain  integration  method  and  the 
DRBEM. 

2.1  Cell  domain  integration 

In  the  cell  domain  technique,  the  attachment  region  Qr  is  subdivided  in  elementary  cells.  The  distri¬ 
bution  of  the  shear  stress  Tj{x.')  in  the  adhesive  is  described  in  terms  of  nodal  values  associated  to  each 
cell.  In  this  work  two  types  of  cells  were  used.  Since  there  exist  two  coincident  nodes  at  crack  elements 
(one  for  each  crack  edge),  these  nodes  can’t  be  used  as  collocation  points  because  no  coincidents  nodes 
exist  in  the  patch.  Then,  constants  cells  with  a  central  node  has  been  used  to  aproximate  the  shear 
stress  distribution  at  neibourghood  of  the  crack.  Nine  node  quadrilateral  isoparametric  cells  were  used 
to  approximate  the  variation  of  the  adhesive  shear  stress  in  the  remaining  attachment  area. 

Consequently,  in  the  cell  integration  method,  the  domain  integral  in  the  equation  (1)  can  be 
expressed  as  (see  Salgado  and  Aliabadi  [13]): 


and  the  integration  is  carried  out  on  each  cell.  Using  equation  (4)  and  the  bi-quadratic  isoparametric 
approximation  proposed  in  this  work,  we  can  write: 


1 


Utfbj  (x)  dQk 


(6) 


where,  N  is  the  matrix  of  bi-quadratic  Lagrange  shape  functions  and  =  |u^,u^|  is  the  vector 
of  nodal  displacements  at  cell  k.  In  this  vector,  refers  to  sheet  displacement  at  Qr  and  refers 
to  repair  displacements.  Similar  expression  can  be  obtained  for  domain  integrals  at  equations  (2)  and 
(3). 


In  this  work  the  integral  on  the  right  hand  side  of  equation  (6)  is  evaluated  using  ten-point  Gaussian 
quadrature.  However,  when  the  source  point  x'  is  placed  within  the  cell,  this  integral  becomes  weakly 
singular  which  will  cause  numerical  error  if  Gaussian  quadrature  is  used  directly.  In  this  case  the 
integrand  in  (6)  can  be  regularized  at  the  singular  point  by  substracting  suitable  singular  term,  which 
may  be  treated  separately  as  follow  (see  Young  and  Rooke  [11]): 


1  -1 

I  U:fNjkdnk  =  I  I  {u:fNjkJ-Xijln{R)j]d^dri 

-1-1 


1  -1 

+XijJ  j  J  ln{R)d^dr] 

-1  -1 


(7) 
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where,  R  =  a/(^  —  +  (^  “  VoY  ■  The  second  integral  on  the  right  hand  side  can  be  evalnated 

analytically.  The  constant  Ajj  is  given  by: 


Aij  — 

where  Gs  is  the  shear  modnlns  of  the  sheet. 


1  (3-^), 

IGtt  Gs 


(8) 


2.2  DRBEM  integration  technique 

In  the  DRBEM,  interaction  forces  are  approximated  as  a  sum  of  unknown  coefficients  af  multiplied 
by  approximating  functions  /j^^(x',x),  so  that: 

D 

bji^)  =  (9) 

The  cofficients  af  have  no  physical  meaning.  But  they  are  related  to  attachment  shear  forces  through 
equation  (4): 


u‘]  (x')  -  uf  (x')  =  J  =  2  (10) 

d=l 

In  this  work,  a  linear  approximation  function  /j^^(x',x)  was  used  for  the  isotropic  sheet: 

//fc  =  (1 (11) 

For  the  anisotropic  patch,  an  approximation  function  given  by  Albuquerque,  Sollero  and  Aliabadi 
[18]  was  used: 

fjk  ~  ^jilm  \cv  (t ,i^lk  T  ^im^lk)]  (12) 

Finally,  the  domain  integral  of  equation  (5)  can  be  expressed  as: 


J  U*f  (x',  x)  bj  (x)  dn^  =  -^J2°‘k  [cij  (x'^)  uij  (x'^)  + 

j  T*f  (x',  x)  ufjdTR  -  j  U*f  (x',  x)  ifjdVR 

Tb  Tr 


(13) 


where  uf^  and  are  particular  solutions  for  displacements  and  tractions  corresponding  to  a  pre¬ 
defined  function  for  the  sheet.  A  similar  approach  was  used  to  model  body  forces  in  the  patch. 


3  Matrix  formulation 


3.1  Cell  integration  technique 

In  matrix  form,  equation  (6)  can  be  written  as: 


1 


\ 

U*  •  NdQk 


k 


a*  -  Ec  Urf  -  Ec  u 


(14) 
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Discretizing  the  boundary,  the  equations  for  isotropic  sheet  (including  traction  equation)  can  be 
written  in  compact  form  as: 

S  _  f^S^S  I  j?S  S  T^S  R 

Me  Ue  -  Me  tc  +  « 

lu^  +  H^uf  =  Gf  tf  +  F^uf  -  (15) 

where  subindex  c  and  d  identify  boundary  and  domain  collocation  points  on  the  sheet.  The  matrix  of 
influence  coefficients  H'^  and  O'®  are  defined  as: 


nelem  p 

=  E  / 

e=l  r 

A  e 

nelem  p 

=  E  y  U*f<l>,dTe  (16) 

e—  1  -p 
A  e 

In  these  integrals,  (f)j  are  shape  functions  for  the  elements.  In  this  work,  quadratic  discontinues 
elements  are  used  to  interpolate  the  displacement  and  traction  variations  in  the  boundaries  of  the 
plate  and  the  repair. 

In  a  similar  way,  matrix  equations  for  repair  can  be  written  as  (without  considering  traction  forces 
applied  at  boundary  repair): 


=  F^^u^  -  F^^u^ 


Iu,^  +  H«  =  F,V-F,V  (17) 

In  this  case,  similar  significance  has  the  and  matrices  as  those  in  the  sheet  case.  In  the 
general  case,  when  the  sheet  and  the  patch  are  made  of  different  materials,  the  F'®  and  F^  matrices 
in  equations  (15)  and  (17)  are  not  equals. 

After  some  mathematical  manipulation,  the  coupling  equations  for  the  sheet  and  the  repair  using 
the  cell  integration  technique  can  be  written  as: 


F'^  1  J  1  _  J  1 

I  J  “  I  0  J 


where  M^,  and  matrix  involving  the  F  matrices  for  sheet  and  repair. 


(18) 


3.2  DRBEM  integration  technique 

In  DRBEM  integration  teenhique,  equation  (13)  can  be  write  in  matrix  form  as: 

y  U*f  (x',  x)  bj  (x)  d^R  =  (19) 

Qr 

In  this  equation,  the  influence  matrices  H'®  and  are  those  defined  in  equation  (16)  with  functions 
ufj  and  t'j^  -  approximated  within  each  boundary  element  by  using  interpolation  functions  and  nodal 
values  as  done  for  (x)  and  tj  (x)  in  equation  (15). 

Discretizing  the  boundary,  equations  for  the  sheet  (including  traction  equation)  can  be  written  in 
a  compact  form  as: 

Hfuf-Gftf  =  Af«^ 


T,,S  I  ’U'S  S  _  A  s  s 

lu^  +  Md  Uc  —  a 


(20) 
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Figure  2:  Model  of  cracked  sheet  reapired  with  adhesive  patch  using  bi-quadratic  interpolation  cells.  Left:  cell 
model.  Right:  DRBEM  model 


where  matrix  A*  is  given  by:  A*  =  In  similar  way,  equations  for  repair  are: 


S 


Hfuf-Gftf  =  Afa^ 

luf +  Hjuf  =  Afa^  (21) 

Now,  equation  (10)  can  be  written  in  a  matrix  form  for  the  sheet  and  the  repair  as: 


ul-n^  = 


-  uf  = 


(22) 


Finally,  coupling  equations  for  the  sheet  and  the  repair  using  the  DRBEM  integration  technique 
are  given  by, 


(h-af-1)^  (af-^)^  1 


(23) 


4  Numerical  results 

4.1  Circular  composite  patch  over  a  cracked  square  sheet 


A  square  sheet  whose  edge  length  is  200  mm  is  subjected  to  a  uniform  constant  tension  of  1  GPa  in 
the  direction  of  the  y-axis.  The  sheet  has  a  central  crack  of  length  2a  =  30  mm  and  thickness  equal 
to  1.5  mm.  A  circular  repair  of  radius  equal  to  30  mm  and  thickness  equal  to  1.5  mm  is  bonded  at 
the  center  of  the  sheet  using  an  adhesive  with  0.15  mm  of  thickness  and  shear  modulus  G  =  0.6  GPa. 
Properties  of  the  sheet  and  the  patch  are  given  in  Table  1. 
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Figure  3:  Normalized  shear  stress  force  in  the  adhesive. 

Table  1:  Mechanical  properties  of  the  sheet  and  the  composite  patch 


Sheet 

Patch 

Young  modulus  {E)  =  72400  Mpa 
Poissons  ratio(t))  =  0.3 

El  =  25000  MPa 

E2  =  208000  MPa 
Gi2  =  72400  MPa 
r'12  =  0.02 

The  problem  was  analyzed  using  the  method  of  cells  and  the  DRBEM.  In  both  cases,  the  mesh 
comprises  of  28  discontinuous  quadratic  elements  on  the  edge  of  the  plate  and  on  the  edge  of  the 
repair.  As  shown  in  Figure  2,  quadratic  continues  cells  with  nine  nodes  were  used  to  discretize  de  load 
transfer  domain  between  the  sheet  and  the  patch  except  in  crack  neighborhood,  where  constants  cells 
were  used.  Ten-point  Gauss  quadrature  rule  was  used  to  evaluate  the  domain  integral  at  quadratic 
cells. 

Also,  figure  2  shows  the  used  DRBEM  model.  In  this  model,  DRBEM  collocations  points  have 
been  concentrated  near  the  crack  and  towards  bonndary  repair.  The  shear  stress  distribntion  in 
the  adhesive  layer  obtained  using  the  DRBEM  is  shown  in  Figure  3.  As  was  expected,  shear  stress 
gradients  appear  near  crack’s  border  where  the  difference  between  sheet  and  repair  displacements  is 
higher.  Shear  distribution  map  obtained  in  the  model  with  cells  is  similar  and  it’s  not  show  here. 

The  resultant  for  the  shear  stress  in  the  adhesive  is  showed  in  the  Figure  4  normalized  with  respect 
to  the  sheet  far  field  stresses  (i.e.  1  GPa).  This  stresses  has  been  obtained  using  the  equation: 

T*  =  ^^\l'^zx  +  Tzy  (24) 

where  (Tq  is  the  far  stresses  applied  in  the  y-axis,  Tzx  and  Tzy  are  shear  stresses  in  the  x  and  y-axis 
directions.  As  can  be  seen  in  this  figure  the  convergence  of  the  solution  is  obtained  as  the  number  of 
internal  points  increases.  Further  refining  in  the  boundary  mesh  hasn’t  significantly  affects  the  results. 
Obtained  results  are  compared  with  analytical  solution  given  by  Rose  [16]  for  an  infinity  orthotropic 
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Figure  4:  Normalized  shear  stress  in  the  adhesive  layer  x=0  and  0  <  y  <  R  <  1. 

patches  bonded  to  an  infinity  orthotropic  sheet  for  patch  with  elliptic  (circular)  geometry: 

T  (y)  =  (25) 

again,  (Tq  is  the  stress  applied  in  the  y-axis  (i.e.  1  GPa)  and  the  parameter  A  is  given  by: 

A2  =  (GA/hA)  I  {E^hs)  '  +  {E^hu)  '  I  (26) 

It  can  be  seen  that  good  agreement  was  obtained  even  for  relatively  coarse  internal  points  grids 
when  the  DRBEM  were  used.  Lower  convergence  rate  to  Rose’s  solution  was  found  with  cell  method. 

4.2  Rectangular  orthotropic  patch  over  a  square  sheet 

Consider  a  thin  aluminium  sheet  with  height  Hg  =  of  254mm,  width  Wg  =  254  mm,  thickness  equal  to 
5  mm  with  a  central  crack  of  length  2a  =  13  mm  repaired  with  boron-epoxi  patch  having  dimensions: 
Wr  =  130  mm;  Hr  =  75  mm.  The  sheet  is  subjected  to  a  remote  uniaxial  tensile  load  of  cr  =  70 
MPa,  plane  stress  condition  are  assumed.  The  material  properties  of  the  plate,  patch  and  adhesive 
are  showed  in  table  2. 

Table  2:  Mechanical  properties  of  the  sheet  and  the  composite  patch 

Sheet  Patch 

Young  modulus  (E)  =  72000  Mpa  Ei  =  19600  MPa 

Poissons  ratio(D)  =  0.33  E2  =  210000  MPa 

Gi2  =  5460  MPa 
1^12  =  0.3 

The  problem  was  analyzed  using  the  cell  method.  The  mesh  comprises  of  28  discontinuous 
quadratic  elements  on  the  edge  of  the  plate.  A  convergence  analysis  for  shear  stress  in  the  adhe¬ 
sive  layer  as  function  of  number  of  cells  and  elements  at  boundary  of  the  repair  was  performed.  Figure 
5  shows  the  used  model. 
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Figure  5:  DEEM  model  for  square  sheet  with  rectangular  patch. 


Figure  6  shows  the  shear  stress  distribution  in  the  adhesive  layer.  Again,  shear  stress  gradients 
appear  near  crack’s  border  where  the  difference  between  sheet  and  repair  displacements  is  higher. 

The  displacement  extrapolation  technique  is  used  for  the  evaluation  of  stress  intensity  factors  as 
described  in  Salgado  and  Aliabadi  [13].  When  discontinues  elements  are  nsed  for  modelling  crack 
surfaces,  SIF  values  are  extrapolated  to  the  crack  tip  using  relationship  (see  figure  7): 

^ -  f (27) 

rAA'  -  rsB'  V  ^aa'  ) 

Three  cases  were  considered,  with  2a  =  13,  15  and  20mm,  respectively.  Table  3  shows  the  stress 
intensity  factors  in  mode  I  obtained  with  12  quadratic  discontinues  boundary  elements  on  each  surface 
of  the  crack.  In  this  table,  SIFs  are  compared  with  those  reported  in  Belhouari  et  al.  [19]. 

Table  3:  KI  stress  intensity  factor  for  rectangular  orthotropic  patch  over  a  square  sheet 


2a(mm) 

Ki-BBM 

A:/-Ref.[19] 

error 

13 

7.60 

8.10 

6.17% 

15 

11.30 

11.90 

5.04% 

20 

11.95 

12.50 

4.40% 

5  Conclusions 

A  new  boundary  element  formnlation  for  modelling  cracked  sheets  repaired  with  composite  patches 
was  developed.  The  cracked  sheet  was  modelled  with  the  DEEM  and  the  patch  was  modelled  with 
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0.15 


0.1 


Figure  6:  Normalized  shear  stress  in  the  adhesive  layer. 


Figure  7:  Discontinuous  crack  tip  element. 


the  BEM.  The  interaction  between  the  isotropic  sheet  and  the  patch  was  modeled  considering  shear 
body  forces  uniformly  distributed  on  the  interaction  zone  using  a  linear  elastic  relationship.  The  cell 
domain  integration  and  the  dual  reciprocity  have  been  used  to  treat  the  domain  integrals  that  arise 
in  the  formnlation  dne  to  shear  interaction  forces.  The  DRBEM  showed  higher  convergence  rate  to 
analytical  solution  than  the  cell  method.  It  can  be  concluded  that  the  new  formulation  can  be  used 
with  reasonable  accuracy  to  study  the  mechanical  behaviour  of  adhesively  bonded  repairs. 
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In  this  work,  an  inverse  problem  of  damage  identification  and  localization  is  modeled  using  independent 
techniques,  both  for  the  direct  model  and  for  the  inverse  model.  The  damage  is  characterized  by  a  hole  in  the 
structure,  which  modifies  existing  temperature  and  stress  fields.  Direct  models  for  the  thermal  (conduction) 
and  elastostatics  problems  are  needed  to  obtain  the  distribution  of  these  quantities  in  the  domain,  for  a 
particular  configuration.  The  boundary  element  method  is  used  here  as  the  direct  problem,  and  two  different 
and  independent  techniques  are  used  for  the  inverse  problem,  in  order  to  localize  and  to  identify  a  damage  in  a 
structure.  The  first  technique  adopted  is  a  global  optimization  technique  using  genetic  algorithms  and  the 
second  approach  is  a  parameter  identification  technique  using  artificial  neural  networks.  The  identification 
and  localization  of  a  hole  in  the  structure  is  performed  using  these  two  techniques  for  the  inverse  problem, 
with  comparable  results. 

Keywords;  damage  detection;  boundary  element  method;  optimization;  genetic  algorithm;  parameter 
identification;  artificial  neural  network 
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1,  Introduction 


Several  types  of  static  and  dynamic  loads  and  the  structural  deterioration  process  can  cause  different 
types  of  structural  damage.  The  damage  can  be  characterized  by  a  change  in  the  structure,  such  as  the 
presence  of  holes  and  cracks.  The  knowledge  of  the  change  in  the  material  properties  corresponding  to  the 
damage  depends  on  the  type  of  material  and  on  the  structural  configuration.  The  proper  assessment  of  the 
damage  in  a  structure  can  be  useful  to  infer  its  remaining  service  life.  The  assessment  of  the  structural  damage 
can  be  performed  through  a  comparison  between  measured  and  simulated  data.  To  provide  the  simulated  data, 
a  numerical  code  is  required,  in  which  a  direct  model  of  the  problem  is  consistently  used  by  an  inverse 
problem  algorithm.  For  the  direct  problem,  a  model  is  required  to  obtain  information  on  the  distribution  of  the 
quantity  of  interest  throughout  the  structure,  given  the  boundary  conditions  and  the  presence  of  the  damage. 
For  the  inverse  problem,  a  model  is  required  for  the  procedure  of  locating  the  damage  in  the  structure  given 
some  (partial)  information  on  the  quantity  of  interest  at  some  particular  locations  (for  example,  where  some 
sensors  are  placed). 

Numerical  methods,  such  as  the  boundary  element  method  (BEM)  or  the  finite  element  method  (FEM) 
can  be  used  for  modeling  the  direct  problem.  The  study  or  analysis  of  damage  in  a  plate  can  be  done  through 
the  thermal  modeling  or  the  distribution  of  stresses.  In  this  work,  two  BEM  formulations  were  used,  for 
potential  and  elastostatics  problems,  respectively.  The  BEM  was  chosen  in  this  work,  for  two  main  reasons;  i) 
the  problem  dimension  under  study  is  reduced  by  one,  and  only  meshing  of  the  boundary  is  required, 
simplifying  the  process  of  re-meshing  the  domain  for  the  various  damages  being  simulated;  ii)  the  integral 
representation  is  an  exact  formulation,  and  the  numerical  errors  are  only  due  to  the  boundary  discretization 
into  boundary  elements.  Eor  the  potential  formulation,  the  potential  values  are  simulated  on  the  external 
surface  of  the  plate  at  given  points.  These  potential  values  represent  the  distribution  of  temperatures  on  the 
plate.  The  use  of  thermal  techniques  shows  that  the  distribution  of  temperatures  on  a  plate  changes  due  to  the 
variations  in  the  mechanical  properties  of  the  plate,  what  could  be  related  to  the  presence  of  a  given  damage. 
Eor  the  elastostatics  formulation,  the  quantities  of  interest  are  the  interior  point  displacements  and  stresses. 

The  damage  detection  problem  can  be  considered  as  a  problem  of  system  identification  or  an  inverse 
problem.  The  inverse  problem  of  identifying  the  presence,  location  and  size  of  damage,  such  as  cracks  and 
holes,  in  a  plate  structure  can  be  modeled  using  optimization  and  parameter  identification  techniques.  In  this 
work,  the  genetic  algorithm  (GA)  is  used  as  the  optimization  procedure,  for  two  reasons:  i)  the  algorithm 
looks  for  a  global  optimum,  and  is  not  trapped  in  local  optima,  which  may  not  locate  properly  the  damage;  ii) 
there  is  no  need  to  evaluate  derivatives  of  the  objective  function.  Also,  in  this  work,  the  artificial  neural 
network  (ANN)  approach  is  used  as  a  parameter  identification  technique,  for  three  reasons;  i)  ANN  does  not 
require  a  priori  the  presence  of  a  Gaussian  white  noise,  as  it  would  be  the  case  of  Kalman  filters,  for  example; 
ii)  ANN  is  capable  of  representing  non-linear  problems;  Hi)  ANN  provides  flexibility  in  terms  of  the  number 
of  internal  layers  to  be  used.  By  solving  the  inverse  problem  using  two  independent  techniques  (GA  and 
ANN),  a  more  reliable  information  on  the  damage  parameters  can  be  obtained,  as  a  comparison  of  the  results 
from  both  approaches  can  provide  a  means  to  verify  these  results. 

The  presence  of  damage  may  induce  rapid  changes  in  the  field  variable  of  the  problem,  and  even 
discontinuities  in  the  governing  equation  in  the  domain.  Classical  calculus-based  optimization  methods 
require  evaluation  of  derivatives  of  the  objective  function,  which  may  not  be  possible  to  be  obtained,  or  may 
be  numerically  obtained,  with  unacceptable  inaccuracy.  Besides,  these  problems  can  have  several  local 
minima  (multiple  solutions),  and  thus  a  global  optimization  method  (such  as  GA)  is  a  better  choice  for  the 
numerical  solution  [1,2].  In  [1],  the  direct  mechanical  problem  is  modeled  numerically  through  BEM,  and  the 
inverse  problem,  to  minimize  the  error  (difference  between  the  measured  and  the  computed  value),  is  modeled 
in  two  ways:  using  sequential  quadratic  programming  (SQP),  to  obtain  a  local  optimum,  and  using  the  GA,  to 
obtain  a  global  optimum  for  the  same  objective  function.  SQP  is  a  calculus-based  optimization  method,  in 


which  the  second  derivatives  are  required  to  obtain  the  Hessian  matrix.  In  [1],  the  Hessian  matrix  was 
approximated  through  a  finite  difference  scheme.  Also,  this  method  depends  on  the  choice  of  the  starting 
point,  so  the  algorithm  can  stop  at  a  local  minimum  of  the  function  that  may  not  represent  the  proper  damage 
parameters  for  the  problem.  On  the  other  hand,  GA  uses  multiple  points  to  search  for  the  solution,  rather  than 
a  single  point,  and  a  global  minimum  has  a  better  chance  of  being  obtained.  Also,  as  GA  does  not  require  any 
evaluation  of  derivatives,  no  errors  are  included  in  the  solution  due  to  the  approximation  of  these  derivatives. 

In  the  works  presented  in  [3]  and  [4],  the  BEM  is  also  used  to  model  the  direct  mechanical  problem 
numerically.  A  backpropagation  neural  network  (BPN)  for  the  on-line  identification  of  holes  or  cracks  in 
composite  structures  is  applied  by  [3].  In  [4],  evolutionary  algorithms  at  the  identification  of  crack  are  used 
and  the  problem  is  formulated  as  the  minimization  of  the  difference  between  the  measured  and  computed 
values  of  displacements  or  stresses  for  selected  boundary  nodes.  The  work  presented  in  [5]  proposes  a  BPN 
for  the  inverse  analysis,  and  a  numerical  model  for  the  direct  method.  This  direct  model  is  based  on  a  coupling 
of  the  FEM  with  the  boundary  integral  equation  (BIE)  method,  of  which  the  discretized  form  is  also  known  as 
the  BEM.  The  BPN  uses  a  backpropagation  learning  rule,  where  the  adjustment  of  weights,  from  input  to 
hidden  layers,  is  made  by  back-propagating  the  errors  of  the  neurons,  from  the  output  layer  to  the  hidden 
layers. 

In  [6],  a  new  method  was  developed  for  finding  boundary  temperatures  and  heat  fluxes,  where  both 
quantities  may  be  unknown  in  some  parts  of  the  boundary.  This  technique  requires  over-specified  thermal 
boundary  conditions,  i.e.,  both  the  temperatures  and  the  heat  fluxes  must  be  specified,  on  other  parts  of  the 
boundary.  In  this  case,  BEM  was  used  for  the  direct  model,  and  the  program  performed  automatic  non¬ 
iterative  determination  of  the  thermal  boundary  conditions  (boundary  temperature  and  flux)  on  the  parts  of  the 
interior  and  exterior  boundaries  where  both  quantities  were  unknown.  This  non-iterative  approach  was 
extended  in  [7]  for  elastostatics  problems  using  BEM,  for  finding  deformations  and  tractions  on  parts  of  the 
boundary  where  these  quantities  are  unavailable.  Again,  the  boundary  conditions  need  to  be  over-specified  on 
other  parts  of  the  boundary,  i.e.,  both  the  displacements  and  the  tractions  must  be  specified  at  these  other 
boundary  subregions. 

In  this  work,  two  direct  problems  are  modeled  using  the  BEM  approach  in  2D;  z)  a  potential  problem 
of  heat  transfer  (conduction)  on  a  domain;  and  ii)  an  elastostatics  problem.  In  both  cases,  a  damage  is 
simulated  by  the  presence  of  a  hole  inside  the  domain.  For  each  run  of  the  direct  model,  the  information  about 
the  location  and  radius  of  the  hole,  and  also  about  the  boundary  conditions,  loading,  and  plate  and  hole 
discretization,  is  also  provided.  After  evaluating  the  boundary  solution,  the  BEM  code  evaluates,  as  a  post¬ 
processing,  some  quantities  of  interest  at  selected  interior  points.  The  selected  interior  points  are  candidates  to 
be  sensor  locations,  for  a  future  experimental  setting,  and  the  quantities  of  interest  at  these  points  may  be  the 
quantities  that  these  sensors  are  able  to  measure.  Each  run  of  the  direct  method  using  the  potential  formulation 
provides  one  piece  of  information  (the  potential,  i.e.,  the  temperature)  at  the  selected  interior  points.  On  the 
other  hand,  the  elastostatics  BEM  formulation  provides  three  pieces  of  information  at  an  interior  point  -  the 
components  of  the  stress  tensor,  i.e.,  two  normal  stresses  and  one  shear  stress.  The  values  of  the  normal  stress 
and  the  shear  stresses  depend  on  the  system  of  coordinates  being  used,  or  on  the  normal  direction  of  the 
cutting  plane  that  passes  through  the  point  of  interest.  As  the  goal  of  the  inverse  method  is  to  identify  and 
locate  the  hole,  but  not  to  identify  any  direction-dependent  properties,  the  desired  quantities  to  be  supplied  to 
the  inverse  model  should  be  scalar  quantities  obtained  at  the  selected  interior  points,  and  not  direction- 
dependent  quantities.  Scalar  quantities  of  interest  can  be  obtained  as  the  invariants  of  the  stress  tensor  -  in  2D, 
the  mean  stress  and  the  octahedral  stress  -  at  the  selected  interior  points.  The  mean  stress  and  the  octahedral 
stress  are  independent  scalar  fields,  and  either  one  can  be  used  as  the  variable  of  interest  at  the  selected 
interior  points.  In  this  work,  the  mean  stress  was  adopted  as  the  quantity  to  be  provided  to  the  inverse  model 
for  the  elastostatics  problem. 


The  boundary  conditions  for  the  external  boundary  of  the  plate  may  be  set  as  temperatures  or  fluxes 
prescribed,  for  the  potential  formulation,  or  displacement  or  traction  prescribed,  for  the  elastostatics 
formulation.  The  boundary  conditions  for  the  internal  boundary  of  the  plate  (the  hole)  were  set  assuming  zero 
fluxes,  for  the  potential  formulation,  and  zero  tractions,  for  the  elastostatics  formulation.  For  the  inverse 
problem,  the  direct  BEM  model  first  evaluates  the  differences  in  the  quantity  of  interest  (the  potential  or  the 
mean  stress,  depending  on  the  problem)  between  the  undamaged  plate  and  the  plate  with  the  damage,  for  all 
selected  interior  points.  These  differences  are  then  supplied  as  input  to  the  optimization  (GA)  or  identification 
(ANN)  subroutines.  The  main  idea  for  passing  only  differences  of  the  quantities  of  interest  is  to  avoid  any 
possible  bias  related  to  the  magnitude  of  these  quantities,  as  only  their  change  (due  to  the  presence  of  the 
hole)  is  important  for  the  inverse  problem.  The  information  provided  by  the  BEM  model  for  the  direct 
problem  is  used  for  comparison  with  similar  information,  which  must  be  available,  for  a  plate  with  a  hole  with 
unknown  size  and  location.  Usually,  the  information  on  the  “real”  plate  would  be  available  by  means  of  an 
experimental  device,  in  which  sensors  would  be  put  in  all  selected  interior  point  locations.  For  the  purpose  of 
validating  this  approach,  the  plate  with  the  “real”  hole  is  also  simulated  with  the  BEM  model,  so  the  inverse 
problem  algorithm  will  try  to  identify  and  locate  this  simulated  “real”  hole.  The  optimization  (GA)  and 
identification  (ANN)  subroutines  are  independent  approaches  for  localization  (obtaining  the  center 
coordinates)  and  identification  (obtaining  the  radius)  of  a  given  simulated  “real”  hole. 

In  short,  the  numerical  modelling  of  the  direct  problem  is  performed  using  two  different  BEM 
approaches,  for  potential  and  elastostatics  formulations,  respectively.  Also,  two  different  and  independent 
techniques  (optimization  using  GA,  and  identification  using  ANN)  are  used  for  resolving  the  inverse  problem 
to  obtain  the  damage  location  and  size,  for  each  BEM  model.  The  comparison  between  the  results  obtained 
using  the  two  different  and  independent  techniques  (GA  and  ANN)  for  the  inverse  problem  allows  for  a 
validation  of  the  inverse  procedure.  A  redundancy  in  the  results,  i.e.,  similar  damage  identification  and 
localization  results,  from  the  two  different  and  independent  inverse  techniques,  will  provide  a  good  indication 
of  the  correctness  of  this  procedure.  A  test  case  using  GA,  available  in  the  literature,  will  also  be  used  for 
comparison  purposes.  All  subroutines  in  this  work  were  written  using  the  MATLAB®  platform. 

2.  Direct  problem:  boundary  element  methods 

Numerical  methods,  such  as  the  boundary  element  method  (BEM)  or  the  finite  element  method  (FEM) 
can  be  used  for  modeling  the  direct  problem.  In  the  FEM,  the  problem  domain  is  partitioned  into  a  number  of 
subdomains  (or  finite  elements)  with  connectivity  between  the  elements  provided  through  common  nodal 
points.  In  the  BEM,  the  governing  partial  differential  equation  of  a  domain  is  transformed  into  a  set  of  integral 
equations,  which  relate  the  boundary  variables  (both  known  and  unknown)  [8,9].  The  BEM  has  some 
advantages  with  regard  to  FEM  [8]:  i)  BEM  discretization  is  done  only  in  the  boundary  of  the  domain,  while 
FEM  requires  the  discretization  of  the  entire  domain;  ii)  the  number  of  equations  associated  with  BEM  is 
smaller  than  in  the  FEM  approach,  for  the  same  degree  of  accuracy;  in)  BEM  is  well  suited  for  problems  with 
singularities,  such  as  in  linear  elastic  fracture  mechanics. 

The  BEM  is  a  numerical  procedure  well  adapted  for  the  modeling  of  a  structure  with  damage.  In  this 
method,  the  distribution  of  the  quantities  of  interest  in  the  domain  is  obtained  from  the  information  of  the 
distribution  of  certain  quantities  in  the  boundary.  Thus,  the  problem  is  described  based  on  what  happens  in  its 
boundaries,  reducing  the  dimension  of  the  problem  and  simplifying  numerically  the  treatment.  In  this  work, 
the  models  investigated  include  the  potential  and  elastostatics  formulations  (see  references  [9]  and  [10]  for 
both  formulations).  A  simple  direct  method  for  a  conduction  problem  is  modeled,  where  the  temperature 
distribution  on  the  external  surface  of  a  thin  plate  is  analyzed.  Without  the  hole,  the  distribution  of  the 
potential  is  known  a  priori.  If  a  small  hole  is  included,  the  potential  distribution  is  unknown  and  must  be 
obtained  numerically  from  the  BEM  solution.  Increasing  the  problem  complexity,  a  BEM  model  for  the 


elastostatics  problem  can  be  used.  Similarly,  the  distribution  of  the  displacement  and  stresses  without  the  hole 
is  known  a  priori.  If  a  small  hole  is  included,  this  information  is  unknown  and  must  be  obtained  numerically 
from  the  BEM  solution.  When  modeling  the  damage  detection  problem  by  means  of  an  analysis  of  the  elastic 
response  of  the  structure  under  excitation,  perturbations  in  the  expected  response  imply  in  the  presence  of 
damage.  Thus,  the  damage  in  the  structure  will  characterize  its  behavior,  static  or  dynamic. 


2.1.  Boundary  integral  equation  for  potential  and  elastostatics  problems 


For  the  elastostatics  problem,  the  elastic  behavior  of  a  body  under  static  loads  is  governed  by  the 
equilibrium,  compatibility  and  constitutive  equations  [11].  Considering  T  as  the  boundary  of  the  body,  an 
integral  representation  of  these  equations,  can  be  written  as  Equation  (1),  for  the  case  with  no  volume  forces 

{y)uk  (j)  =  \X^‘k  (^)  -  (x)}/r(x)  (i) 

where;  q^.  is  the  traction  vector  at  a  boundary  point  whose  outward  normal  n . ;  is  the  displacement  vector; 


u[  and  q[  are  the  displacement  and  traction  vectors  of  the  fundamental  solution,  respectively  (see  references, 
[11],  [7],  [9]  and  [10]  for  more  details).  When  the  limit  to  the  boundary  is  taken  for  the  collocation  point  y , 
the  equation  is  called  a  Boundary  Integral  Equation  (BIE).  The  term  c[  is  the  free  term  coefficient,  which 
depends  on  the  position  of  the  collocation  point,  relative  to  the  boundary.  For  an  interior  point,  c[  =  1 ;  for  an 
exterior  point,  =  0;  for  a  boundary  point  on  a  smooth  section  of  the  boundary,  c[  -  Xjl .  For  non-smooth 
boundary  points,  the  free  term  coefficient  depends  on  the  swept  angle  at  this  point,  when  going  from  the 
boundary  region  before  the  point  to  the  boundary  region  after  the  point,  following  the  interior  domain.  By 
performing  collocation  at  different  boundary  points  (the  nodes),  a  set  of  equations  is  obtained,  which  can  be 
discretized  to  obtain  a  system  of  algebraic  equations  to  be  solved.  The  set  of  equations  is  completed  by  the 
boundary  conditions,  ufx)^u.  on  r„  and  ^,(x)  =  ^,.  on  T^,  where  T,,  and  T^  are  non-overlapping  partitions 

of  the  boundary  T  (T^uT^^T  and  r„nr^=0)  ([11,12]).  The  kernels  of  the  integrands,  given  by  the 


fundamental  solution  and  its  derivative,  lead  to  weakly-singular  and  singular  integrals,  respectively,  when  the 
collocation  point  and  the  integration  point  coincide.  Special  integration  schemes  are  incorporated  in  the  BEM 
code,  to  account  for  the  evaluation  of  these  singular  integrals. 

Equation  (1)  is  a  component  of  a  vector  equation,  in  the  k  -direction  (k  =  1,2  ,  in  the  2D  case).  A  scalar 
boundary  integral  equation  for  the  potential  problem  can  be  obtained  as  an  integral  representation,  closely 
similar  to  Equation  (1),  for  the  case  where  there  are  no  heat  sources  in  the  domain.  In  this  case,  due  to  the 
scalar  nature  of  the  potential  field,  the  symbol  k  can  be  dropped,  as  shown  in  Equation  (2). 


c{y)u{y)^\^ 


^(x;j)^(x) - [x;y^u[x^  c/r(x) 

dn 


(2) 


where:  c[y)  is  the  coefficient  of  the  free  term;  u  is  the  potential;  q  is  the  flux  in  the  outward  normal 
direction;  ^  =  (l/2;r)ln(l/r)  is  the  fundamental  solution  for  the  Eaplace  equation;  r-\x-y\  is  the  distance 
between  the  collocation  point  y  and  the  integration  point  x;  and,  dy//dn--{\l2;i:r)[drldn)  is  the  flux 
associated  to  the  potential  y/  .  The  boundary  conditions  are  similar  to  the  previous  case,  with  u  and  q  now 
representing  known  values  of  the  potential  and  flux  on  r„  and  T^ ,  respectively  [12]. 


2.2.  Boundary  element  discretization 


By  evaluating  Equation  (1)  at  the  collocation  points  y,  by  using  proper  shape  functions  in  the 


discretized  boundary  (in  this  work,  eonstant  boundary  elements),  and  by  applying  adequate  quadrature 
formulae  for  the  numerical  integration  (in  this  work,  4  Gauss  points  for  eaeh  element),  a  system  of  linear 
equation  is  obtained  as  (Equation  (3)) 

[H\{u]^[G]{q]  (3) 

where  [u]  and  [q]  contains  the  nodal  values  of  the  displacement  and  traction  vectors,  for  the  elastostatics 
problem,  or  the  nodal  values  of  potential  and  flux  vectors,  for  the  potential  problem. 

When  the  boundary  eonditions  of  eaeh  problem  are  taken  into  aeeount  properly,  after  algebraic 
manipulation,  known  and  unknown  quantities  are  separated,  and  a  system  of  linear  equations,  whieh  can  be 
solved  for  the  unknown  boundary  quantities,  is  obtained  as  (Equation  (4)) 

[A]{x]^{f]  (4) 

where  {x}  is  the  vector  of  unknown  boundary  quantities;  {/}  is  the  right-hand  side,  obtained  after 
manipulating  the  known  boundary  quantities  with  the  proper  numerieal  integration  eoefficients;  and  [A\  is  a 
matrix  with  the  integration  eoefficients  related  to  the  unknown  boundary  variables. 

After  the  boundary  solution  is  obtained,  by  post-processing,  the  solution  for  the  displacement 
(elastostaties  problem)  and  for  the  temperature  (potential  problem)  at  seleeted  interior  points  is  obtained  by 
means  of  a  partieular  ease  of  Equation  (1),  where  c[  is  equal  to  1  [10].  As  the  integral  equation  for  interior 
points  does  not  eontain  singular  integrals,  special  integration  sehemes  are  not  required  in  the  BEM  eode,  for 
this  ease. 

Regarding  the  elastostaties  problem,  the  internal  stresses  ean  be  eomputed  by  differentiating  the 


displacements  at  internal  points  and  introducing  the  eorresponding  strains  into  the  stress-strain  relationships 
(Equation  (5))  (see  references  [9]  and  [10]). 
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where;  //  is  the  shear  modulus  and  v  is  the  Poisson’s  ratio. 

After  a  proper  substitution  of  the  value  of  into  Equation  (5),  the  internal  stresses  ean  be  represented  in 
eompaet  form  shown  in  Equation  (6) 
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In  Equation  (7),  «  =  1 ,  ^2,  and  y  =  4  for  a  two-dimensional  ease.  The  derivatives  indicated  by  commas  are 

taken  at  a  boundary  point  xf  (Equation  (8)) 
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with  rf  =  xf  -  x/  ,  being  x/  an  internal  point  and  =  |r| .  The  derivative  shown  in  Equation  (8)  is  equal  and 
opposite  in  sign  to  those  taken  at  an  internal  point  [9]. 


3,  Inverse  problem:  optimization  and  parameter  identification  techniques 

The  inverse  problem  might  be  modeled  by  means  of  optimization  and  parameter  identification 
techniques.  The  damage  is  simulated  by  the  presence  of  a  small  hole  in  the  domain,  and  the  goal  is  to  obtain 
size  and  location  of  the  damage. 

3.1.  Optimization  using  genetic  algorithms 

The  genetic  algorithm  (GA)  is  a  search  method  based  on  the  processes  of  natural  evolution.  This 
method  works  with  a  set  of  possible  solutions  for  a  given  problem,  composing  the  initial  population.  In  other 
words,  GA  uses  multiple  points  to  search  for  the  solution  rather  than  a  single  point  in  the  traditional  gradient 
based  optimization  method  [13].  In  this  algorithm  the  problem  variables  are  represented  as  genes  in  a 
chromosome  (each  chromosome  is  also  denominated  an  individual  of  the  population).  Starting  from  an  initial 
population,  the  individuals  with  better  adapted  genetic  characteristics  have  higher  chances  of  surviving  and 
reproducing. 

According  to  [4],  the  GA's  are  methods  that  do  not  depend  on  the  choice  of  the  initial  point,  increasing 
the  chances  of  obtaining  the  optimum  global  of  the  system.  So  that  the  population  is  diversified  and  maintain 
certain  acquired  adaptation  characteristics  by  the  previous  generations,  the  genetic  operators  (selection, 
crossover  and  mutation)  can  be  used.  These  operators  transform  the  population  through  successive 
generations,  extending  the  search  until  arriving  to  a  satisfactory  result.  For  more  details  about  how  these 
operators  work,  see  references  [14],  [15]  and  [16]. 

In  this  work,  the  optimal  solution  for  unknown  parameters  of  the  damage  (location  and  size)  is 
obtained  through  the  GA  for  potential  and  elastostatics  formulation.  Considering  the  first  formulation 
(potential  formulation),  a  functional  can  be  defined  as  the  difference  between  the  measured  (simulated)  values 
of  the  difference  in  the  potential  (between  the  undamaged  plate  and  the  plate  with  the  damage)  and  the  values 
of  the  same  differences  in  potential  calculated  at  the  same  points  by  the  damage  detection  program.  In  the 
second  formulation  (elastostatics  formulation),  the  functional  is  defined  as  the  difference  between  the 
measured  (simulated)  values  of  the  local  difference  in  the  mean  stress  (between  the  undamaged  plate  and  the 
plate  with  the  damage)  and  the  values  of  the  same  differences  in  mean  stress  calculated  at  the  same  points  by 
the  code  (assuming  several  different  locations  and  sizes  for  the  “numerical”  damage).  The  functional 
corresponds  to  the  fitness  function  of  the  GA.  The  minimization  of  this  fitness  function  allows  the  damage 
detection  program  to  find  the  unknown  parameters  of  the  damage.  The  functional  formulation  is  shown  at 
Equation  (9): 

1  ” 
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being  n  the  number  of  internal  points  i  (“sensors”  placed  in  the  plate)  where  the  differences  are  evaluated; 
measuredi  the  vector  of  simulated  values  for  the  differences  obtained  using  BEM,  for  a  given  damage;  and, 
calculatedji  the  vector  of  differences  in  potential  (potential  formulation)  or  mean  stress  (elastostatics 
formulation)  calculated  by  the  code  for  each  individual  j. 

As  mentioned,  GA  starts  with  an  initial  population,  representing  a  set  of  possible  solutions  for  a  given 
problem.  To  solve  the  damage  detection  problem,  each  chromosome  (individual)  of  the  population  can  be 
assembled  according  to  the  vector  presented  in  Equation  (10); 

C  =  g2  ^3  (10) 

where; 

gi  -  first  gene  representing  the  x-coordinate  of  the  hole; 
g2  -  second  gene  representing  the  y-coordinate  of  the  hole; 


g3  -  third  gene  representing  the  hole  radius; 

g4  •••  gn+3  -  fourth  gene  and  the  subsequent  genes,  representing  the  measures  of  the  potential  difference 
(potential  formulation)  or  the  mean  stress  difference  (elastostatics  formulation)  between  the  undamaged  plate 
and  the  plate  with  the  damage. 

As  an  example,  Figure  1  represents  three  possible  configurations  of  chromosomes.  While  the  location 
and  size  of  the  hole  vary,  the  number  and  location  of  the  sensors  remain  the  same,  for  all  chromosomes.  The 
information  on  the  quantity  of  interest  is  collected  at  these  sensor  locations  for  all  cases. 


Figure  1.  Plate  with  a  hole:  three  possible  eonfigurations  for  the  ehromosomes. 

3.2,  Parameter  identification  using  artificial  neural  network 

The  artificial  neural  networks  (ANN's)  are  computational  techniques  that  present  a  mathematic  model 
to  represent  the  human  brain  and  to  try  to  simulate  the  learning  process  of  this  brain.  An  ANN  is  formed  by 
the  interconnected  neurons  whose  inputs  can  be  obtained  from  the  outputs  of  other  neurons  or  from  input 
nodes.  Different  configurations  of  the  artificial  neuron  can  be  made  to  develop  different  network  topologies 
[18].  Among  the  existent  configurations,  the  ANN  can  be  feedforward  or  feedback.  At  the  feedforward  neural 
networks,  the  neurons  are  interconnected  in  layers,  but  the  flow  of  data  only  occurs  in  a  direction  [17].  At  the 
feedback  neural  networks,  a  neuron  receives  the  information  of  neurons  of  the  previous  layer  and  of  a 
subsequent  layer.  After  defining  the  structure  of  the  ANN,  an  iterative  process  of  weight  adjustment  of  this 
network  is  made.  This  process  is  known  as  training  process.  Following  the  training,  the  ANN  learns  how  to 
proceed  for  other  input  data  in  the  problem  domain. 

In  this  work  a  backpropagation  neural  network  (BPN)  is  used,  through  a  feedforward  configuration 
and  the  backpropagation  learning  algorithm.  The  backpropagation  algorithm  carries  out  a  supervised  learning 
where  the  desired  outputs  are  given  as  part  of  the  training  vector.  For  more  details  about  how  this  algorithm 
woks,  see  reference  [19].  In  addition,  there  is  a  training  function  of  ANN  known  as  ‘gradient  descent  with 
momentum  and  adaptive  learning  rate  function’  (see,  for  example  the  MATLAB®  help  manual).  This  function 
is  a  backpropagation  network  training  function  that  combines  adaptive  learning  rate  with  momentum  training. 
An  adaptive  learning  rate  allows  the  performance  of  the  steepest  descent  algorithm  to  improve,  attempting  to 
keep  the  learning  step  size  as  large  as  possible  while  keeping  learning  stable.  Moreover,  momentum  training 
allows  a  network  to  respond  not  only  to  the  local  gradient,  but  also  to  recent  trends  in  the  error  surface. 
Without  momentum  a  network  may  get  stuck  in  a  flat  local  minimum. 

3.3.  Formulation  of  optimization  and  parameter  identification  problems 

The  problem  of  damage  detection  in  a  thin  plate  can  be  formularized  as  an  optimization  problem 
(using  GA)  according  to  the  flowchart  in  Figure  2,  or  can  be  formularized  as  a  parameter  identification 
problem  (using  ANN)  according  to  the  flowchart  in  Figure  3. 

Considering  Figure  2,  the  initial  population  for  the  GA  approach  is  formed  by  the  geometric 
information  of  a  numerical  hole  (x-  and  y-coordinates  of  its  center,  and  also  its  radius)  and  also  by  differences 
in  the  quantities  of  interest,  calculated  at  selected  interior  points,  herein  called  “Difference  1”.  “Difference  1” 
can  be  the  local  difference  in  the  potential  or  the  local  difference  in  the  mean  stress  between  the  undamaged 


plate  and  the  plate  with  the  damage,  for  potential  and  elastostaties  formulations,  respectively.  Similarly,  a  set 
called  “Difference  2”  can  be  evaluated  at  the  same  interior  points,  representing  the  “measured”  differences  for 
the  quantity  of  interest  at  these  points,  for  the  “real”  hole.  In  this  work,  the  “real”  hole  is  also  simulated.  To 
validate  the  damage  detection  approach,  the  value  of  “Difference  2”  was  not  allowed  to  be  in  the  initial 
population  of  the  GA  approach.  The  initial  population  and  also  “Difference  2”  are  employed  in  the  fitness 
function,  presented  in  Equation  (9).  The  goal  of  the  GA  approach  is  to  look  for  a  minimum  value  of  this 
fitness  function.  For  that,  the  algorithm  uses  genetic  operators  to  modify  the  population  and  subsequently 
reevaluate  the  fitness  function  for  the  new  population.  Convergence  criteria  can  be  set,  including  the  number 
of  iterations,  the  differences  in  the  fitness  function  (the  optimum  value)  between  two  subsequent  populations, 
or  the  differences  between  the  hole  parameters  of  location  and  size  (the  optimizer)  between  two  populations. 
When  the  convergence  criterion  is  met,  the  numerical  holes  have  reached  the  vicinity  of  the  “real”  hole,  and 
thus  the  information  about  the  location  and  size  of  the  “real”  hole  is  obtained.  In  this  work,  a  criterion  for  the 
maximum  number  of  generations  (no  to  exceed  75  in  the  potential  problem  and  not  to  exceed  100  in  the 
elastostaties  problem)  was  assumed,  together  with  a  default  criterion  for  the  tolerance  (difference  between  two 
fitness  functions  less  than  or  equal  to  1  x  10“*’ ). 


Figure  2.  Flowchart  for  the  optimization  procedure  using  GA. 

According  to  Figure  3,  a  network  is  created,  considering  “Difference  1”  (the  same  “Difference  1”  as  in 
the  GA  approach)  as  the  input  data  and  the  geometric  information  for  the  hole  (x-  and  y-coordinates  of  the 
hole  center,  and  its  radius)  as  the  output  data.  The  next  step  is  to  train  the  created  network,  obtaining,  as  a 
result,  a  NET  that  contains  information  about  how  to  proceed  for  another  input  data  in  the  problem  domain. 
Finally,  the  trained  network  is  simulated  for  “Difference  2”  (same  “Difference  2”  as  in  the  GA  approach). 
Similarly  to  the  optimization  algorithm,  convergence  criteria  need  to  be  set  for  this  approach.  In  this  work,  the 
error  goal  was  assumed,  not  to  exceed  1x10“^  for  the  potential  problem  and  1x10“^  for  the  elastostaties 


problem.  The  order  of  magnitude  of  these  assumed  error  goals  follows  the  order  of  magnitudes  in  the 
differenees  of  the  quantities  of  interest,  namely  the  potential  and  the  stress.  A  eonvergenee  eriterion  in  terms 
of  the  maximum  number  of  iterations  (not  to  exceed  5000  epochs)  was  also  assumed.  When  the  convergence 
criterion  is  met,  the  ANN  has  identified  the  “real”  hole,  providing  the  information  about  its  location  and  size. 


Figure  3.  Flowchart  for  the  parameter  identification  procedure  using  ANN. 

4.  Numerical  results  and  discussion 


For  the  potential  problem,  the  results  obtained  by  the  damage  detection  program  are  analyzed  for  a 
problem  of  heat  flow  in  a  thin  plate.  Initially,  a  plate  without  damage  and  with  the  dimensions  (0.06x0.06)  m 

was  simulated  through  the  boundary  element  method  (BEM),  as  illustrated  in  Figure  4.  The  boundary  of  the 
plate  was  discretized  into  12  elements  and  the  value  of  the  potential  was  evaluated  at  49  internal  points 
(Figure  4(b)).  The  contour  conditions  for  the  problem  are  represented  in  this  figure,  where  q  represents  the 
heat  flow  and  u  represents  the  temperature  at  the  boundary.  Then,  a  plate  with  a  central  hole  of  radius 
0.06  cm  ,  with  the  same  dimensions  and  boundary  conditions,  was  also  simulated,  and  the  obtained  results  for 
the  potential  were  compared  with  the  plate  without  damage. 


</  =  o 
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Figure  4.  Plate  model  for  the  potential  problem:  (a)  dimensions,  loading,  and  boundary  conditions,  (b)  boundary  discretization  and 

sensor  locations. 


For  the  elastostatics  problem,  a  BEM  model  was  built  for  the  plate  with  a  hole  with  the  boundary 
eonditions  illustrated  in  Figure  5(a).  Two  diseretizations  were  implemented  for  the  external  contour,  a  coarse 
mesh  with  12  constant  elements  and  a  fine  mesh  with  48  constant  elements.  Figure  5(b)  shows  the 
discretization  for  the  case  of  48  elements  in  the  outer  boundary  and  12  elements  in  the  hole,  as  well  as  the 
position  of  the  nine  sensors.  At  the  present  work,  the  sensors  were  uniformly  distributed  on  the  plate  and  no 
positioning  study  of  the  sensors  was  performed.  The  plate  was  simulated  with  shear  modulus  equal  to 
94,500  MPa  and  a  Poisson’s  ratio  for  plane  strain  equal  to  0. 1 . 


Figure  5.  Plate  model  for  the  elastostatics  problem:  (a)  dimensions,  loading,  and  boundary  conditions.  Insert  shows  a  stress-free 
hole;  (b)  boundary  discretization  (fine  mesh)  and  sensor  locations.  Insert  shows  hole  discretization. 

The  influence  of  the  numerical  errors  due  to  the  BEM  discretization  in  the  optimization  results  can  be 
seen  in  Eigure  6.  A  comparison  was  made  for  the  optimization  results  (using  10  runs  of  a  GA  approach)  and 
for  two  meshes  (a  coarse  mesh  and  a  fine  mesh)  using  the  elastostatics  formulation  for  the  plate  shown  in 
Figure  5(a).  Figure  6  presents  illustrative  results  for  the  mean  values  of  the  error  in  the  location  (x  and  y 
coordinates)  and  size  (radius  r)  of  a  central  hole  [20]. 


Figure  6.  Mean  values  for  the  error  in  location  and  size  of  a  center  hole. 

4.1.  Assembly  of  data  for  the  GA  and  ANN  procedures 


For  the  potential  formulation,  holes  with  radius  equal  to  0.15  cm,  0.03  cm  and  0.09  cm,  respectively. 


were  considered  to  assemble  the  initial  population  of  the  GA.  For  each  one  of  these  radius,  the  coordinate  x 
of  the  center  of  the  hole  was  varied  from  0.5  cm  to  5.0  cm  and  the  coordinate  y  of  the  center  of  the  same  hole 
was  varied  from  0.5  cm  to  5.5  cm,  both  coordinates  with  a  step  size  of  0.5  cm.  Then,  110  different  positions 
for  each  radius  in  the  plate  were  simulated  and  the  respective  values  of  the  potential  difference  were  stored  for 
a-posteriori  processing. 

For  the  ANN,  initially  25  internal  points,  representing  the  sensors  on  the  plate,  were  considered  to  the 
assembly  of  the  input  data  of  this  network.  After,  the  number  of  sensors  was  subsequently  decreased  to  15,  9 
and  5,  respectively.  In  the  present  work,  the  sensors  were  uniformly  distributed  on  the  plate  and  no 
positioning  study  of  the  sensors  was  performed;  only  a  study  regarding  the  reduction  on  number  of  the  sensors 
was  done.  The  distribution  of  the  sensors  on  the  plate,  for  each  case,  is  shown  in  Figures  7(a),  7(b),  7(c)  and 
7(d),  respectively. 


(c) 
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Figure  7.  Sensor  distribution  for  the  ANN:  (a)  25,  (b)  15,  (c)  9  and  (d)  5. 


In  this  potential  formulation,  a  single  hole  with  a  radius  equal  to  0.15  cm  in  nine  different  hole 
positions  was  considered  to  assemble  the  input  and  output  data  of  ANN.  (positions  equal  to  (0.5;0.5)  cm, 
(0.5;3)  cm,  (0.5;5.5)  cm,  (3;0.5)  cm,  (3;3)  cm,  (3;5.5)  cm,  (5.5;0.5)  cm,  (5.5;3)  cm  and  (5.5;5.5)  cm).  Then, 
another  hole  of  radius  equal  to  0.05  cm  was  also  analyzed  in  each  mentioned  position. 

For  the  elastostatics  formulation,  holes  with  radius  equal  to  0.05  cm,  0.10  cm  and  0.15  cm  was 
considered  to  assemble  the  data  GA  and  ANN.  For  each  radius,  the  x  and  y-coordinate  of  the  center  of  the 
hole  was  varied  from  0.5  cm  to  5.5  cm  with  a  step  size  of  0.5  cm.  Then,  121  different  positions  for  each  radius 
in  the  plate  were  simulated  and  the  respective  values  of  the  difference  in  the  mean  stress  at  the  9  internal 
points  (shown  as  sensor  locations  in  Figure  5(b))  were  found  by  means  of  BEM  and  these  values  were  stored 
for  a-posteriori  processing. 

4.2.  Analysis  of  the  results  obtained  from  the  genetic  algorithm 

For  the  potential  formulation,  in  the  initial  population  of  the  GA,  the  values  of  the  difference  in  the 


potential  were  normalized  before  using  the  data  directly,  taking  into  consideration  the  maximum  value  of  this 
difference.  Finally,  the  initial  population  with  330  individuals  can  be  formed.  As  the  potential  values  near  the 
right  border  (temperature  equal  to  zero)  of  the  plate  are  close  to  zero,  the  potential  difference  is  used  instead 
of  the  direct  use  of  the  potential  value. 

The  plots  of  the  location  and  size  of  the  holes  obtained  from  5  different  runs  of  the  GA  are  presented 
in  Figure  8.  The  program  was  run  only  5  times,  because  there  was  no  significant  difference  when  this  value 
was  increased.  The  “real”  position  of  the  hole  is  represented  in  continuous  line  and  the  results  found  by  the 
GA  in  non  continuous  lines. 


Figure  8.  “Real”  and  simulated  hole  for  potential:  (a)  for  a  eentral  hole  with  elitism  equal  to  2;  (b)  for  a  eentral  hole  with  elitism 
equal  to  10;  (e)  for  a  hole  at  (4;5)  cm.  Inserts  shows  the  region  of  hole  in  detail. 

As  presented  previously  (in  section  3.3),  a  stopping  criterion  was  assumed  for  the  number  of 
generations  no  to  exceed  75  in  the  potential  problem,  and  a  stopping  criterion  for  the  tolerance  (difference 
between  two  fitness  functions  less  than  or  equal  to  1x10“'’).  In  the  results  obtained,  the  tolerance  of  the 
problem  was  reached,  in  other  words,  there  was  no  improvement  in  the  objective  function  (fitness  function), 
and  the  maximum  number  of  generations  was  not  reached,  showing  a  good  convergence  of  the  algorithm.  The 
crossover  fraction  was  set  as  0.8,  and  the  mutation  fraction  was  set  as  0.2.  Figure  8  (a)  and  Figure  8  (b)  show 
a  hole  in  the  position  (3;3)  cm,  and  Figure  8  (c)  shows  a  hole  in  the  position  (4;5)  cm,  all  cases  considering  a 
radius  equal  to  0.06  cm.  For  the  first  simulation  (Figure  8  (a)),  the  elitism  was  2,  and  the  second  (Figure  8  (b)) 
and  third  simulation  (Figure  8  (c)),  the  elitism  was  changed  from  2  to  10,  guaranteeing  that  10  individuals 
survive  in  the  next  generation.  With  the  change  of  elitism  (Figure  8  (b))  the  holes  were  concentric,  and  the 
hole  position  presented  a  small  uncertainty.  Moreover,  the  radius  for  every  simulation  was  not  much  sensitive 
to  the  variation  of  the  GA  parameters.  The  crossover  function  considered  was  the  heuristic  function  with  a 
value  of  ratio  equal  to  1.3  (this  value  represents  how  far  the  child  is  from  the  better  parent).  Besides,  the 
mutation  function  adopted  was  the  Gaussian  function.  The  results  are  different  for  each  run  of  GA  approach, 
because  there  is  a  small  mutation  presence  and  a  crossover  function  that  is  different  for  each  run  of  the 
algorithm,  in  other  words,  there  is  an  associated  occurrence  probability. 

In  a  similar  manner,  for  the  elastostatics  formulation,  the  values  of  the  difference  in  the  mean  stress 
were  normalized,  taking  into  consideration  the  maximum  value  of  this  difference.  The  values  of  x  and  y- 
coordinate  of  the  center  of  the  hole  and  its  radius  were  also  normalized,  considering  the  respective  maximum 
values.  After  that,  the  initial  population  with  363  individuals  can  be  formed. 

The  plots  of  the  location  and  size  of  the  holes  obtained  from  10  different  runs  of  the  GA  are  presented 
in  Figure  9.  The  GA,  due  to  its  own  randomness,  generates  a  different  optimal  solution  every  time  it  is  run; 
nevertheless  the  results  of  the  GA  approach  present  a  tendency  to  be  concentrated  near  the  “real”  hole.  Figure 


9  (a)  shows  the  results  for  a  eentral  hole;  Figure  9  (b)  shows  a  hole  located  at  (2;2)  cm;  and  Figure  9  (c),  a 
hole  located  at  (5;3)  cm.  The  radius  of  each  plot  was  considered  equal  to  0,12  cm.  The  “real”  position  of  the 
hole  is  represented  in  continuous  line  and  the  results  found  by  the  GA  in  non  continuous  lines. 


(a)  (b)  (c) 

Figure  9.  “Real”  and  simulated  hole  for  mean  stress:  (a)  for  a  central  hole;  (b)  for  a  hole  at  (2;2)  cm;  (c)  for  a  hole  at  (5;3)  cm. 

Stopping  criteria  were  assumed  in  the  elastostatics  problem,  both  for  the  number  of  generations  (no  to 
exceed  100),  and  also  for  the  tolerance  (difference  between  two  fitness  functions  less  than  or  equal  to  1x10“'’). 
The  crossover  fraction  was  set  as  0.95,  and,  hence,  the  mutation  fraction  is  0.05  for  the  GA  approach 
presented  in  this  work.  The  elitism  considered  that  10  individuals  survive  in  the  next  generation.  The  function 
that  performs  the  crossover  was  heuristic  function,  considering  a  value  of  ratio  equal  to  0.9.  The  mutation 
function  was  uniform  function  where  each  gene  has  a  probability  0.03  of  being  mutated.  Another  GA 
parameter  configured  was  the  migration.  In  this  work,  the  migration  fraction  was  set  as  0.20,  the  direction  of 
migration  was  set  as  “both”  directions  and  20  generations  pass  between  migrations  of  individuals  between 
subpopulations.  Finally,  the  selection  function  was  the  roulette  selection. 

The  GA  technique  requires  soma  extra  care  for  its  implementation,  due  to  the  required  choices  for  the 
configuration  of  the  algorithm  parameters,  which  may  be  different  for  each  problem.  This  choice  depends  on 
the  realization  of  a  great  number  of  experiments  and  tests.  Moreover,  the  GA  also  presents  a  high 
computational  cost  due  to  the  several  evaluations  of  the  fitness  function.  The  damage  detection  code  using  GA 
can  find  a  region  for  the  probable  occurrence  of  the  hole,  as  this  algorithm  generates  a  different  optimal 
solution  every  time  it  is  run.  Thus,  a  confidence  interval,  for  the  different  parameters  being  identified,  can  be 
obtained. 

4.3.  Analysis  of  the  results  obtained  from  the  artificial  neural  network 

Considering  the  problem  of  heat  flow,  initially  the  presence  of  a  single  hole  in  the  structure  was  studied. 
Then,  the  influence  in  the  results  was  verified  when  the  number  of  sensors  at  the  plate  was  decreased. 
However,  as  above  mentioned,  no  study  regarding  the  sensor  positioning  was  accomplished  in  this  work.  As 
well  as  the  input  data,  the  values  used  to  test  the  network  were  assembled  following  the  sensor  distribution 
scheme  at  the  plate.  A  hole  of  radius  0.10  cm  in  different  positions  was  considered  to  test  the  network.  The 
best  choice  for  the  parameters  of  the  backpropagation  neural  network  (BPN)  was  50  neurons  in  the  input 
layer,  4  neurons  in  the  hidden  layer,  and  4  neurons  in  the  output  layer.  The  other  parameters  of  the  ANN  were 
set  as: 

•  Threshold  function  in  the  input  and  hidden  layers;  tan-sigmoid  transfer  function; 

•  Threshold  function  in  the  output  layer;  linear  transfer  function; 

•  Training  function;  gradient  descent  with  momentum  and  adaptive  learning  rate; 


Error  goal:  1x10“^; 
Number  of  epochs:  5000; 
Learning  rate:  0.05. 


The  influence  of  the  reduction  of  the  sensor  number  in  the  results  found  by  the  ANN  for  a  hole  in  the 
position  (3;3)  cm  and  radius  equal  to  0.10  cm  can  be  analyzed  at  Table  1.  Table  2  shown  the  results  for  a  hole 
in  the  position  (4;2)  cm.  The  problem  domain  is  reduced  when  there  is  a  decrease  of  the  sensor  number  on  the 
plate.  The  obtained  results  depend  on  the  distribution  of  the  sensor  on  the  plate  and  of  the  quality  of  the  input 
data. 


Table  1.  Influence  of  the  reduction  of  the  sensor  number  for  a  central  hole. 


“Real”  hole 

Simulated  hole 

Sensors  number 

X 

y 

r 

X 

y 

r 

25 

3.00 

3.00 

0.10 

3.0035 

3.0003 

0.0992 

15 

3.00 

3.00 

0.10 

2.9949 

2.9977 

0.1009 

9 

3.00 

3.00 

0.10 

2.9998 

2.9973 

0.1002 

5 

3.00 

3.00 

0.10 

3.0010 

2.9959 

0.1000 

Table  2.  Influence  of  the  reduction  of  the  sensor  number  for  a  non-central  hole. 


“Real”  hole 

Simulated  hole 

Sensors  number 

X 

y 

r 

X 

y 

r 

25 

4.00 

2.00 

0.10 

3,4568 

0,5676 

0,0994 

15 

4.00 

2.00 

0.10 

2,1225 

0,5135 

0,0994 

9 

4.00 

2.00 

0.10 

2,4224 

0,4355 

0,0224 

5 

4.00 

2.00 

0.10 

1,4138 

0,9774 

0,1000 

The  results  obtained  for  a  hole  of  radius  0.10  cm  in  the  positions  (3;3)  cm  (Figure  10(a)),  (1;1)  cm 
(Figure  10  (b)),  and  (5;5)  cm  (Figure  10  (c))  for  5  sensors  on  the  plate  is  shown  in  Figure  10. 


(a)  (b)  (c) 

Figure  10.  Potential  problem:  results  from  the  ANN  with  5  sensors  for  a  hole  at  position:  (a)  (3;3)  cm;  (b)  (1;1)  cm;  (c)  hole  at 

position  (5;5)  cm. 

For  the  elastostatics  formulation,  the  ANN  simulates  the  non-linear  behavior  between  the  values  of  the 
local  difference  in  the  mean  stress  (between  the  undamaged  plate  and  the  plate  with  the  damage)  and  the  hole 
parameters  (location  and  size).  Information  regarding  the  difference  in  the  mean  stress  is  supplied  in  the  input 
of  the  network,  besides  the  parameters  of  the  hole  are  supplied  in  the  output  of  the  same  network.  Holes  of 


different  sizes  and  at  different  places  can  be  part  of  the  data  supplied  to  the  net.  Having  defined  the  input  and 
output  data,  the  next  step  is  to  build  the  network  and,  then,  this  network  can  be  trained.  Finally,  the  network 
can  be  tested  for  other  data  of  difference  in  the  mean  stress,  obtaining  as  answer,  the  location  and  size  of  the 
hole. 

As  according  to  the  initial  population  of  GA,  the  values  of  the  difference  in  the  mean  stress  and  the 
hole  parameters  were  also  normalized,  before  using  these  values  directly.  After  training  the  network  with 
these  data,  this  network  was  tested  for  a  hole  of  radius  0.12  cm  in  different  positions.  Figure  11  shows  some 
the  results  obtained,  considering  9  sensors  on  the  plate  and,  whose  distribution  is  presented  in  Figure  5(b). 
The  network  was  configured  with  100  neurons  in  the  input  layer,  50  neurons  in  the  hidden  layer,  and  3 
neurons  in  the  output  layer.  Differently  from  the  parameter  configuration  of  ANN  for  potential  formulation, 
the  error  goal  was  set  as  1  x  10“\  and  the  learning  rate  was  set  as  0.01 . 


(b) 


(c) 


Figure  11.  Elastostatics  problem:  results  from  the  ANN  for  a  hole  at  position:  (a)  (3;3)  cm;  (b)  (1;1)  cm;  (c)  (5;5)  cm;  and. 


In  Figurell,  the  results  present  a  small  area  of  uncertainty  near  the  “real”  hole;  moreover,  the  size  of 
the  hole  was  obtained  with  good  accuracy.  These  results  were  similar  to  those  presented  in  Figure  10,  for  the 
potential  problem,  and  they  were  obtained  more  quickly  than  in  the  case  of  using  GA  (as  a  global  optimization 
technique).  For  this  reason,  the  solution  of  a  damage  detection  problem  through  the  ANN  (as  a  parameter 
identification  technique)  is  also  known  as  an  online  identification.  An  advantage  of  the  use  of  ANN  in  regard 
to  the  GA  is  that,  after  training  the  network,  holes  with  different  sizes  and  in  different  locations  can  be  tested 
without  running  the  damage  detection  program  again. 

The  damage  detection  problem  using  parameter  identification  technique  was  solved  more  quickly  than 
in  the  case  of  using  global  optimization  techniques.  In  this  work,  the  solution  of  the  problem  through  ANN 
presented  good  results  for  the  several  parameters  being  identified.  In  particular,  the  size  of  the  hole  was 
obtained  with  good  accuracy,  and  the  location  of  the  hole  was  given  by  a  fairly  small  area  of  uncertainty  near 
the  “real”  hole,  for  the  several  cases  tested.  In  part,  difficulties  in  finding  the  exact  area  of  the  occurrence  of 
the  damage  are  due  to  training  problems  of  the  network,  or  choice  of  the  configuration  parameters  of  the 
network  or  the  choice  of  the  input  and  output  data.  Taking  into  account  the  advantages  of  each  technique,  a 
hybrid  approach  could  be  considered  for  future  work.  In  this  approach,  the  GA  could  be  used  to  find  the 
occurrence  area  of  the  damage,  and  then  the  ANN  could  find  the  exact  size  of  this  damage,  reducing  the 
search  time  for  the  optimum  result. 


4.4.  Example  of  analysis  of  noise  or  measurement  error  in  the  data 

To  examine  how  the  inverse  method  using  GA  herein  responds  to  measurement  error,  random  noises 
were  introduced  into  measured  data.  The  flowchart  presented  in  Figure  12  shows  this  approach. 

The  random  noise  is  a  signal  formed  by  a  set  of  random  numbers  drawn  from  a  normal  distribution 


with  zero  mean  (white  noise)  and  with  COV  (eoefficient  of  variation)  given  as  a  pereentage  (5%  or  10%)  of 
the  measurement  value  at  the  sensor  location.  This  noise  is  added  to  the  measured  data,  to  create  a  set  called 
“Measured  data  2”.  This  new  measured  data  was  normalized  (as  discussed  in  section  4.2)  and  then  used  in  the 
GA  approach  for  the  elastostatics  problem.  The  GA  approach  was  run  10  times,  for  each  case  (5%  and  10% 
noise),  always  considering  the  same  configuration  of  parameters  as  in  the  case  without  noise.  In  each  run  of 
the  GA,  a  different  noise  signal  was  generated,  with  the  proper  COV. 


Figure  12.  Flowchart  for  the  analysis  of  the  measurement  error. 


A  hole  in  (3,3)  cm  position  with  a  radius  size  equal  to  0.12  cm  was  simulated  for  the  elastostatics 
problem,  considering  a  random  noise  with  5%  and  10%  introduced  into  measured  data.  The  results  are 
summarized  in  Table  3  and  Table  4.  Table  3  shows  the  mean  results  obtained  through  10  runs  of  GA  with  5% 
random  noise  in  measured  data,  where  errors  around  2.10%,  3.30%  and  4.17%  were  obtained  in  x-location,  y- 
location  and  radius  size,  respectively,  which  are  comparable  with  the  errors  for  the  noise-free  results.  These 
errors  were  small  and  no  regularization  is  needed  in  this  case.  Besides,  as  mentioned  previously,  GA  generates 
a  different  optimal  solution  every  time  it  is  run.  When  the  GA  was  analyzed  with  10%  random  noise  (Table 
4),  the  error  in  radius  size  was  about  7.50%,  and  the  errors  in  x-location  and  y-location  were  similar  to  the 
errors  for  the  5%  random  noise  case.  These  results  show  that  the  GA  optimization  procedure,  for 
identification  and  localization  of  the  hole  in  the  structure,  presents  very  small  sensitivities  to  changes  in  the 
measured  values  at  the  sensors,  proving  the  robustness  of  the  algorithm.  A  similar  analysis  of  the 
measurement  noise,  or  errors  in  the  measured  data,  can  be  performed,  to  investigate  the  case  of  the  inverse 
method  using  ANN,  and  is  the  object  of  the  current  research. 


Table  3.  Results  obtained  with  GA  for  5%  random  noises  into  measured  data. 


Result 

noise-free 

Error  (noise-free)  [%] 

With  noise  of  5% 

Error  (with  noise)  [%] 

X 

3.069 

2.30 

3.063 

2.10 

y 

2.731 

8.97 

2.901 

3.30 

r 

0.127 

5.83 

0.125 

4.17 

Table  4.  Results  obtained  with  GA  for  10%  random  noises  into  measured  data. 


Result 

noise-free 

Error  (noise-free)  [%1 

With  noise  of  10% 

Error  (with  noise)  [  %] 

X 

3.069 

2.30 

3.006 

0.20 

y 

2.731 

8.97 

2.933 

2.23 

r 

0.127 

5.83 

0.129 

7.50 

If  a  logarithmic  transformation  j'.  -  log(  ^ )  (with  ^  =  0. 1  x  10  ^  to  prevent  the  appearance  of  a  -go 

value  in  the  function  [1])  is  done  in  the  funetional  formulation  (Equation  (4)),  the  result  for  a  central  hole  is 
slightly  improved.  Table  5  shows  the  comparison  between  the  results  without  logarithmic  transformation 
(case  A)  and  with  this  transformation  (ease  B)  for  a  hole  in  (3,3)  cm  position  and  a  radius  size  equal  to  0.12 
cm,  and  errors  of  simulation  with  regard  to  a  “real”  hole.  Table  6  shows  the  mean  results  obtained  through  10 
runs  of  GA  with  10%  random  noise  in  measured  data,  where  an  error  about  0.57%,  2.23%  and  7.50%  were 
obtained  in  x-location,  y-loeation  and  radius  size,  respeetively,  which  are  eomparable  with  the  errors  for  the 
noise-free  results,  considering  a  logarithmie  transformation. 


Table  5.  Comparison  of  GA  identification  results  using  the  mean  stress,  with  and  without  a  logarithmic  transformation  (cases  B  and 

A,  respectively). 


Result 

Case  A 

Error  case  A  [%] 

Case  B 

Error  case  B  [%] 

A 

3.069 

2.30 

2.918 

2.73 

Y 

2.731 

8.97 

2.909 

3.03 

R 

0.127 

5.83 

0.122 

1.67 

Table  6.  Comparison  of  GA  results  for  the  mean  stress  with  logarithmic  transformation  (case  B),  with  and  without  10%  random 

noise  in  the  measured  data. 


Result 

Case  B 

Error  case  B  [%] 

With  uoise  of  10% 

Error  (with  uoise)  [%1 

X 

2.918 

2.73 

3.017 

0.57 

y 

2.909 

3.03 

2.933 

2.23 

r 

0.122 

1.67 

0.129 

7.50 

The  results  from  Table  5  illustrate  the  fact  that,  in  most  cases,  the  use  of  the  logarithmic 
transformation  tends  to  reduee  the  pereent  error  in  the  identification  of  the  parameters  of  the  damage. 
Aecording  to  Table  5,  the  error  iny-location  for  both  cases,  A  and  B,  was  larger  than  the  error  in  x-location.  A 
possible  reason  for  the  differenee  in  these  results  is  that  only  a  central  hole  was  simulated  for  the  plate  model 
presented  in  Figure  5.  From  Table  6,  one  can  see  that  no  significant  change  oeeurred  in  the  results  when  a 
noise  of  10%  was  added  in  the  measured  data.  The  GA  approaeh  presented  in  this  work  is  robust  in  regard  to 
measurement  error,  as  only  small  errors  were  obtained  at  the  results  (radius,  x-  and  y-loeation)  when  an  error 
of  10%  was  added  in  the  measured  data.  The  plate  is  square  and  symmetrie,  as  ean  be  seen  in  Figure  5; 
however,  the  boundary  conditions  induce  an  asymmetry  in  the  model.  The  influence  of  the  plate's  aspect  ratio 
and  boundary  conditions,  as  well  as  the  proximity  between  the  hole  and  different  sections  of  the  boundary  (for 
example,  smooth  parts  of  the  boundary  versus  corners,  or  sections  with  different  boundary  conditions)  is  an 
object  of  current  investigation. 

The  present  work  is  limited  to  a  simple  thin  plate  with  a  cireular  hole.  For  more  general  problem, 
multiple  damages  ean  be  addressed,  requiring  a  new  assembly  of  ehromosomes  (Equation  (10))  of  GA  and 
also  changes  in  the  input  data  of  ANN,  to  account  for  the  presence  of  these  multiple  damages.  Cracks  can  also 
be  modeled,  on  a  first  approximation,  as  elliptical  holes,  so  that  the  BEM  formulation  presented  in  this  work 
could  also  be  used.  By  doing  so,  a  new  assembly  of  chromosomes  of  the  GA  has  to  be  performed,  including 
individuals  that  represent  elliptical  hole,  and  new  inputs  have  to  given  to  ANN  aecordingly,  allowing  the 
identifieations  of  eircular  and  elliptieal  holes  in  the  plate.  After  this  first  approximation  for  the  craek  is 
implemented,  a  more  aeeurate  formulation  for  the  BEM  direct  model  should  be  considered,  to  properly  model 
the  presence  of  the  erack  (for  example,  by  using  a  BEM  formulation  speeific  for  fraeture  mechanics,  such  as 
the  dual  boundary  element  method  [21]).  As  part  of  the  ongoing  researeh,  the  extension  of  this  approach  to 
multiple  damages  and  multiple  types  of  damages  is  being  undertaken,  together  with  a  proper  treatment  of 


uncertainties,  whieh  are  present  not  only  in  the  measurements,  but  also  in  the  numerieal  simulations  (due  to 
diseretization  errors),  and  in  the  problem  parameters,  sueh  as  the  domain  geometry  variables  and  material 
properties. 

4.5.  GA  approach  for  damage  identification:  a  comparison  with  literature  results 

The  method  herein  (that  uses  measurements  of  differenees  in  mean  stress)  was  eompared  with  a  result 
presented  in  [1]  (that  uses  measurements  of  the  boundary  displaeements  and  traetions).  For  both  examples  that 
were  eompared,  a  plate  with  external  dimensions  (O.lOxO.lO)  m  was  simulated,  and  the  loading  was  applied 
on  the  left-hand  side  external  boundary.  The  material  eonstants  were  eonsidered  equal  to  100  GPa  for  shear 
modulus  and  0.3  for  Poisson’s  ratio.  In  our  work,  the  results  of  eomparison  were  reaehed  for  a  statie  loading 
of  1000  MPa  in  both  the  horizontal  and  in  the  vertical  coordinate  direetion,  on  the  left-hand  side  external 
boundary,  and  the  right-hand  side  was  fixed.  In  [1],  the  plate  was  subjeeted  to  a  harmonie  dynamie  loading  in 
both  direetions  at  the  left-hand  side,  and  the  right-hand  side  was  also  fixed. 

Table  7  shows  the  results  for  GA  found  by  [1]  and  the  results  presents  in  our  work.  In  [I],  the  results 
were  obtained  after  running  GA  with  200  generations,  a  population  equal  to  5  (no  information  is  given  in  that 
text  on  how  the  individuals  of  the  population  are  plaeed  in  the  plate),  and  the  plate  under  dynamie  loading  (for 
more  details,  see  the  referenee  [1]).  In  our  work,  the  parameters  of  the  GA  were  eonfigured  aeeording  to  those 
parameters  presented  in  seetion  4.2  for  the  elastostaties  problem,  regarding  200  generations  and  a  population 
of  49  individuals.  As  a  test,  only  a  hole  with  diameter  equal  to  0.5  was  eonsidered  in  some  positions  where  the 
test  ease  (“real”  hole)  was  not  ineluded  in  the  initial  population  (the  x  and  y-eoordinate  of  the  eenter  of  the 
hole  was  varied  from  0.5  em  to  9.5  em  with  a  step  size  of  1.5  em),  henee,  validating  the  results  obtained.  For  a 
general  problem,  more  individuals  have  to  eonsider  in  the  initial  population  of  GA. 

Table  7.  GA  approach:  comparison  with  literature  results. 


Test 

“Real” 

hole 

Results  presented  by  [1] 

Results  in  this  work 

Calculated  best 
element 

Average  for  1000 
solutions 

Error 

[%1 

Calculated  best 
element 

Average  for  20 
solutions 

Error 

[%1 

X 

4.0 

3.9606 

5.59 

38.75 

3.7336 

3.52 

12.00 

y 

4.0 

4.0236 

4.74 

18.50 

3.9578 

3.95 

1.25 

diameter 

0.5 

0.4968 

0.52 

4.00 

0.5000 

0.53 

6.00 

As  shown  in  Table  7,  the  GA  approaeh  used  in  this  work  has  presented,  for  most  oases,  more  aoourate 
results  in  the  identifioation  of  the  “real”  hole  dimensions,  with  respeot  to  the  GA  approaeh  used  in  [1].  In  the 
literature  example,  an  average  of  1000  solutions  was  oomputed,  while,  in  this  work,  only  an  average  for  20 
solutions  was  performed.  Also,  for  eaoh  solution,  only  a  few  seoonds  were  needed  to  run  the  inverse  program 
using  GA  on  a  PC.  These  features  illustrate  the  aeeuracy  and  the  low  eomputational  eost  of  the  eurrent 
approaeh. 

5.  Conclusions 

In  this  work,  an  inverse  problem  of  identifying  damage  in  a  plate  strueture  was  solved  using  both 
optimization  and  parameter  identifieation  teehniques.  A  genetie  algorithm  (GA)  was  used  as  the  optimization 
teehnique,  and  an  artifieial  neural  network  (ANN)  eode  was  used  as  the  parameter  identifieation  proeedure. 
Two  models  for  the  direet  problem  were  investigated,  one  eonsidering  a  heat  flow  problem  and  another 
eonsidering  an  elastostaties  problem.  In  the  heat  flow  problem,  the  boundary  element  method  (BEM)  for  the 


potential  was  considered  for  the  direct  problem.  The  BEM  for  the  potential  supplies  the  necessary  information 
(potential  values  at  internal  points  of  the  plate)  to  the  damage  detection  program.  In  the  elastostatics  problem, 
a  boundary  element  method  (BEM)  formulation  was  used  as  the  direct  model  in  this  inverse  problem.  The 
refinement  of  the  mesh  for  the  direct  BEM  model  was  shown  to  play  an  important  role,  improving  accuracy  in 
the  damage  identification  results,  when  a  fine  mesh  was  used.  The  analyses  of  the  results  indicates  that  the 
damage  detection  code  using  GA  can  only  find  a  region  for  the  probable  occurrence  of  the  hole,  as  this 
algorithm  generates  a  different  optimal  solution  every  time  it  is  run.  The  fitness  function  of  the  GA  approach 
presented  in  this  work  has  converged  for  the  specified  tolerance,  before  the  algorithm  has  reached  the 
maximum  number  of  generations.  Moreover,  this  GA  approach  was  robust  in  regard  to  the  measurement  error, 
as  only  a  small  error  was  obtained  in  the  results  when  a  noise  of  10%  was  added  to  the  measured  data.  Also, 
this  GA  approach  compares  well,  both  in  accuracy  and  in  computational  cost,  with  respect  to  a  similar  GA 
approach  used  in  the  literature  for  damage  identification.  The  solution  of  the  problem  through  ANN  has  also 
presented  good  results  for  the  several  parameters  being  identified. 

An  important  observation  is  that  very  small  holes  are  difficult  to  observe  by  the  damage  detection 
program,  mainly  when  these  holes  are  close  to  the  borders  of  the  plate.  The  optimization  and  the  identification 
techniques  adopted  in  this  inverse  problem  can  be  used  concomitantly,  as  independent  procedures  to  identify 
the  presence  of  a  hole  on  the  plate,  thus  providing  a  means  to  verify  the  numerical  results  obtained  for  the 
location  and  size  of  the  damage  in  the  structure,  increasing  the  confidence  in  the  damage  identification  results. 
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Abstract 

This  paper  presents  a  boundary  integral  formulation  for  the  computation  of  moments  and  stresses  at  internal  and 
boundary  points  of  laminated  composite  plates.  An  integral  equation  for  the  second  displacement  derivative  is  develo¬ 
ped  and  all  derivatives  of  the  fundamental  solution  are  computed  analytically.  Stresses  on  the  boundary  are  computed 
by  a  procedure  that  uses  integral  equations,  derivatives  of  shape  functions,  and  constitutive  relations.  The  obtained 
results  are  in  good  agreement  with  finite  element  results  available  in  literature. 

Key  words:  laminated  composite,  boundary  element  method,  stress  analysis. 


1.  Introduction 

The  attempt  of  developing  analytical  models  for  the  rep¬ 
resentation  of  the  behavior  of  plates  comes  since  middle 
of  1800  with  works  developed  by  Sophie  Germain,  La¬ 
grange,  and  Poisson  [1].  Since  1978,  when  the  first  gen¬ 
eral  direct  formulation  based  on  the  Kirchhoff’s  hypo¬ 
thesis  appeared,  the  boundary  element  method  (BEM) 
has  had  large  growth,  being  nowadays  applied  to  se¬ 
veral  practical  engineering  problems.  The  first  works 
discussing  the  use  of  boundary  element  direct  formula¬ 
tion,  in  conjunction  with  the  Kirchhoff’s  theory,  were 
by  Bezine  [2],  Stern  [3],  and  Tottenhan  [4].  Nowa¬ 
days,  BEM  is  a  well-established  numerical  technique  to 
deal  with  an  enormous  number  of  engineering  complex 
problems.  Analysis  of  plate  bending  problems  using 
the  BEM  has  attracted  the  attention  of  many  researchers 
during  the  past  years,  proving  to  be  a  particularly  ade¬ 
quate  field  of  applications  for  that  technique.  The  fun¬ 
damental  solution  is  an  essencial  part  of  the  boundary 
element  method.  Bending  analysis  of  thin  plates  by  the 
BEM  requires  the  use  of  two  fundamental  solutions:  the 
displacement  field  due  to  a  transverse  point  load,  and  the 
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displacement  field  due  to  a  point  moment.  Fundamental 
solutions  for  anisotropic  plates  utilize  complex  variable 
theory  following  Lekhnitskii  [5].  Shi  and  Bezine  [6], 
presented  a  boundary  element  analysis  of  plate  bending 
problems  using  fundamental  solutions  proposed  by  [7] 
based  on  Kirchhoff  plate  bending  assumptions.  Rajamo- 
han  and  Raamachandran  [8],  proposed  a  formulation 
where  the  singularities  were  avoided  by  placing  source 
points  outside  the  domain.  Paiva  et  al  [9],  presented  an 
analytical  treatment  for  singular  and  hypersingular  inte¬ 
grals  of  the  formulation  presented  in  [6].  Albuquerque 
et  al  [10]  presented  a  method  to  transform  domain  in¬ 
tegrals  into  boundary  integrals  in  the  formulation  pre¬ 
sented  in  [6].  In  [11],  this  formulation  was  extended 
for  dynamic  problems.  Shear  deformable  shells  have 
been  analyzed  using  the  boundary  element  method  by 
[13]  with  the  analytical  fundamental  solution  proposed 
by  [14].  Wang  and  Huang  [15],  presented  a  boundary 
element  formulation  for  orthotropic  shear  deformable 
plates.  Later,  in  [16],  the  previous  formulation  was 
extended  to  laminate  composite  plates.  Recently,  [17] 
presented  a  displacement  discontinuity  formulation  for 
modeling  cracks  in  orthotropic  Reissner  plates. 

This  paper  proposes  numerical  procedures  to  com¬ 
pute  moments  and  stresses  in  internal  points  and  at  the 
boundary  of  composite  laminated  plates  using  a  boun¬ 
dary  element  thin  plate  formulation.  To  the  best  of  au¬ 
thor’s  knowledge,  there  is  no  paper  in  literature  that 
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presents  a  boundary  element  formulation  to  compute 
moments  and  stresses  in  anisotropic  plates  or  shells. 

2.  Boundary  integral  equation 


the  corners;  P  is  the  held  point;  Q  is  the  source  point; 
and  an  asterisk  denotes  a  fundamental  solution. 

The  transverse  displacement  fundamental  solution  is 
given  by: 


As  shown  by  [10],  boundary  integral  equations  for  j 

transverse  displacements  of  anisotropic  thin  plates  and  ~  ^  {C\R\{p,S)  +  C2R2(p,0) 

its  derivative  can  be  written,  respectively,  as: 

Kw(Q)  +C3  [5i(p,  0)  -  52(p,  0)]} , 
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where  ^  is  the  derivative  in  the  direction  of  the  outward 
vector  n  that  is  normal  to  the  boundary  T;  ^  is  the  di¬ 
rectional  derivative  at  the  source  point  on  the  boundary 
normal  direction  m;  M„  and  y„  are,  respectively,  the 
normal  bending  moment  and  the  Kirchhoff  equivalent 
shear  force  on  the  boundary  T;  7?^  is  the  thin-plate  reac¬ 
tion  of  the  corners;  Wc  is  the  transverse  displacement  of 


+  >  RciiP)^ 
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X  and  y  are  the  coordinates  of  the  held  point  P,  Xg  and 
yo  are  the  coordinates  of  source  point  Q, 
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Figure  1 :  Staking  sequence  in  a  laminate  composite  plate  with  N  plies. 


where  matrix  Q  is  given  by: 

Q  =  T  ‘Q(T 


(15) 


The  transformation  matrix  T  is  given  by: 
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where  a  is  the  angle  between  the  fiber  direction  of  given 
ply  and  the  global  reference  xy,  and  stiffness  matrix  Q 
is  given  in  terms  of  engineering  constants  by: 
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di  and  e,  are  the  real  and  the  imaginary  part,  respec¬ 
tively,  of  the  roots  of  the  characteristic  polynomial: 

D22P^  +  ^D2(,p^  +  2(Di2  +  2D^(,)p^ 


+4Di6jU  +  Dii  =0,  (13) 

Dll,  Di^,  Di2,  D22,  D26,  and  Dgs  are  the  bending  stiff¬ 
ness  constants  of  an  anisotropic  thin  plate. 

The  repeated  index  i  in  /?,  and  S ,  does  not  imply  sum¬ 
mation.  The  coefficient  a  is  an  arbitrary  constant  taken 
as  fl  =  1. 

Other  fundamental  solutions  are  shown  in  works  of  [6] 
and  [10]. 


3.  Computation  of  stresses  and  moments  on  internal 
points 

Laminates  are  fabricated  such  that  they  act  as  an  integral 
structural  element.  To  assure  this  condition,  the  bond 
between  two  plies  in  a  laminate  should  be  infinitesi¬ 
mally  thin  and  not  shear  deformable  to  avoid  the  ply  slip 
over  each  other,  and  to  allow  displacement  continuity 
along  the  bond  (see  [21]).  Thus,  we  could  consider  that 
strains  are  continuous  along  the  thickness.  However,  as 
each  ply  is  of  a  different  material,  stresses  present  dis¬ 
continuities  along  laminate  interfaces. 

As  presented  by  [21],  stresses  at  each  ply  can  be  evalu¬ 
ated  from  strain  by: 
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where  and  Ej  are  elasticity  moduli  in  the  longitudi¬ 
nal  and  transversal  directions,  respectively,  Glt  is  the 
shear  modulus  in  the  plane  of  the  ply,  and  vlt  is  the 
Poisson  coefficient  (see,  for  example  [21]). 

Moments  are  given  by: 
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In  Kirchhoff  plates,  strains  are  given  by: 
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where  z  is  the  distance  from  the  point  where  displace¬ 
ments  are  been  computed  and  the  midplane  (see  Figure 
1). 

So,  in  order  to  obtain  strains,  moments,  and  stresses,  the 
second  order  derivatives  of  integral  equation  (1)  need  to 
be  computed.  These  derivatives  are  given  by: 
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where  the  second  derivative  of  the  transversal  displace¬ 
ment  fundamental  solutions  in  relation  to  x  is  given  by: 
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Second  order  derivatives  of  other  fundamental  solutions 
in  relation  to  x  are  given  by: 


dy^ 

X 


dx^ 


^  d^w*  ^  d^ 
Jl^TX  +/2 


w  .  (9‘*v 

+  /3- 


(9x4  (9x^5y  dx^dy^ 


d'^V*  d^M*  dw(P) 

-^(Q,P)w(P)  -  -^(Q,P)-^ 


dy^ 


dT(P) 


^  d^Rl 

+  Z  ^<2- 

i=l 


5  -R*, 

clx^ 


=  -\81 


d-^w* 


d^w* 


+  82 


(9x4  ^^dx^dy  ^^dx^dy^ 


d^w* 


X 


+  y/^c,(^)^(G-^)+  (25) 

4^  5y2  5y2 


52 


w  ,  „  „  . .  (9^v 


Vn(P)-^(Q,P)-M„(P) 

52 


5n5y' 


(2,7’) 


5x2 


w  ,  5^ 

+  /l2 


w  .  5^ 


,  (28) 


,  (29) 


„  w  5^w  ' 

^3  ^  ^  ,  +/J4: 


5x® 


5x4  5y  5x2  5y  2  dyp-dy^ 


dY(P) 


1  I  d^w*  ,  5V*  ,  5V* 


(30) 


52w(2) 

5x5y 

-I 


52  y*  d^M*  dw(P) 

'’(Q,P)w(P)-^(Q,P) 


dxdy 


dxdy 


dn 


dY(P) 


Second  order  derivatives  in  relation  to  y  and  xy 
are  computed  by  similar  procedures  and  will  not  be 
shown  herein.  Derivatives  of  fundamental  solutions  of 
transversal  displacement  w  can  be  expressed  by  linear 
combination  of  /?,  and  S ,  derivatives.  All  derivatives  of 
Rj  and  5,  up  to  the  4th  order  are  presented  by  [6].  The 
5th  order  derivatives  are  given  by: 
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(31) 


"idi  [d^  -  {dj^  +  sin^  6  cos  6 


Qi j{.  8 [di  cos^  6  +  3  [dj^  +  e,-^)  sin  0  cos^  0 

dx^dy  Ai 

3di  [dj^  +  sin^  6 cos  0  +  [di^  -  sin^  6] 


A2 


{d,^  -  {d,^  +  sin^  0 


,  (32) 


A2 


d^S  i  Aei  sin6[-3  cos^  0  -  6di  sin  6  cos  0 


(36) 


1 

00 

[  [c/,7  -  e,7^  cos^  6  +  3c/;  [c/,7  +  e,7^  sin  6cos^  6 

Ai 

dx^dy^ 

[  A2 

[e/^  - 

1  ^ 

■  3c//^)  sin^  6] 

3  {di^  +  e,7)' 

sin^  6  cos  6  +  c/;  [c/,7  +  e,7^  sin^  6 1 

^  (33) 

Ai 

+ 

A2  j 

d^Si 

4e;[cos^  6-3  [c/,7  +  e,7^  sin^  6  cos  6 

dx^dy 

Ai 

d^Ri  _  8  ^ 

[  di  [c/,7  -  3e,7^  cos^  6  +  3  [c/,7  -  e,7 

)  sin6cos^  6 

dx^dy^  p^ 

1  A2 

2c/;  [c! 

f,7  +  e,7)  sin^  6] 

(37) 


Ai 


(38) 


3di  {d^  +  e,7)  sin^  0  cos  0  +  {d^  +  e,7)  sin^  0 


A2 


’  d^S  i  4e,  [2t/;  cos^  6  +  3  +  e,7)  sin  6  cos^  6 


dx^dy^ 


Ai 


5^/?;  _  8  1 

[c/,7  -  6e,7c/,7  +  e,7j 

1  cos^  6 

1 

[c/,7  +  e,7) 

dxdy"^  p^  1 

A2  -- 

Ai 

(39) 


3di  (^di*  -  2ep'di^  -  36,“*)  sin  6cos^  6 


A2 


,  4ei  cos  0[{3di^  -  e,7)  cos^  6 


3  (di^  -  e,7)  (di^  +  e,7)  sin^  6  cos  6 


di  (di^  +  sin^  6 

A2  I 


(35) 


dx^dy^  Ai 

3  (di^  +  e,7)  sin  6  [2c/;  cos  6  +  [c/,7  +  e,7^  sin  6]  j 
Ai 

d^Si  4e;  f  4c/;(c/;  -  e;)(c/;  +  g;)  COS^  6 


,  (40) 


dxdy* 


A2 


1 

00 

1  Cl' 

';  [c/,7  -  10e,7c/,7  +  5e,7)cos^  6 

3  [3c/,7  -  e,7)  [c/,7  +  e,-2^ 

1 

w  1 

^ ^ 

1 

A2 

^  A2 

3  (c/;®  -  Sei^di^  -  5e,4c/,7  +  ef'^  sin  6cos^  6 
A2 


6c/;  [c/,7  +  e,7^  sin^  6  cos  6  +  [c/,7  +  e,7^  sin^  6 

A2  I 


(41) 
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Figure  2:  Variables  of  the  transformation  of  the  domain  integral  into  Figure  3:  Normal  m  and  tangential  t  directions  at  the  source  point, 

a  boundaiy  integral. 


and 


d^Si  Aci  [  ~  lOe p-di^  +  cos^  6 

dy^  p3  I  ^2 

12di  (di^  -  e,'*)  sin  0cos^  6 
A2 


3  -  3di^^  {di^  +  sin^  0cos  0 

A2 


4.  Computation  of  stresses  on  the  boundary 

Equations  (24),  (25),  and  (26)  present  integrals  whose 
kernels  have  singularities  of  order  To  compute 
stresses  on  the  boundary,  those  equations  will  provide 
integrals  that  are  more  than  hypersingular.  Conse¬ 
quently,  the  computation  of  stresses  on  the  boundary 
demands  an  alternative  approach. 

The  directional  derivative  of  equation  (1)  at  the  source 
point  in  the  boundary  tangential  direction  (Figure  3)  is 
given  by; 


dM^dw 
dt  dn 


dt 


2di{di^  +  sinks'] 


A2 


(42) 


where  Ai  -  [cos^  6  +  {di^  +  sin^  6  +  di  sin  26^ 


l^cos^  6  +  [di^  +  sin^  6  +  di  sin  26j  . 

.r 

Jq 

f  I  dw*  d^WX  ■y-~' 


'’"'■jo 


(45) 


The  last  term  of  equation  (24)  can  be  transformed  from 
a  domain  integral  to  a  boundary  integral  following  the 
procedure  proposed  by  [10].  Thus,  considering  a  linear 
distributed  load  g  -  Ax  +  By  +  C,  we  have; 


Inside  a  quadratic  discontinuous  boundary  element,  the 
directional  derivative  of  the  transversal  displacement  in 
the  tangential  direction  is  given  by; 


d^w’ 
^  dx^ 


dO.  = 


X 


H*  dr 
r  on 


where 


H* 


(Ax  +  By +  C) 


d^w*  , 


(43) 


(44) 


dw 

dt 


Nv 


dw\ 

dt 


+  N2 


dW2 

dt 


-t  A^3 


5w3 

dt 


(46) 


where  and  ^  are  derivatives  of  w  at  boundary 

element  nodes  in  the  tangent  to  the  boundary  direction 
at  the  source  point.  A^i,  N2,  and  Nj,  are  quadratic  dis¬ 
continuous  shape  functions  written  as; 


Ni 


r  is  the  value  of  p  at  the  boundary  (Figure  2). 
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(47) 


Ni 


(48) 


where  D.j  (i,j  =  1,2,6)  are  the  flexural  rigidities  in 
reference  system  mt,  that  is: 


and 


A^3 


(49) 


D'n 

Di2 

D  = 

^12 

^22 

^26 

II 

H 

6 

H 

(57) 

i  ^'l6 

D'2, 

D',e  \ 

The  second  derivative  of  the  transversal  displacement  is 
given  by: 

where 

Du 

D\2 

Die 

d^w  d  ldw\d^ 

D  = 

D\2 

D22 

D26 

'dfi  ~  ^  ^ 

D\6 

D26 

Dee 

(58) 


Writing  ^  in  terms  of  nodal  values  interpolated  by 
shape  functions  as  given  by  equation  (46),  we  have: 


In  this  system,  the  known  variables  are:  M„„  ^  and 

^2  f<2 

and  the  unknown  variables  are:  M,,  M^t,  and 


d^w 

dfi 


dw 


dw 


5w3  \  d^ 


(51) 


dt 


dt 


dt  I  dt 


d^w  I  dN\  dw\  dN2  dw2  dN^  dw^  ) 
dfi  \  d^  dt  d^  dt  d^  dt  j 


(52) 


Following  a  similar  procedure  to  it  is  possible  to 
obtain: 


d^w  I  dNi  dw\  dN2  dw2  dN^  dw^ 
dmdt  \  d^  dm  d^  dm  d^  dm 


X 


(53) 


where  and  ^  are  derivatives  of  w  at  boundary 

element  nodes  in  the  normal  to  the  boundary  direction 
at  the  source  point. 

It  was  not  possible  to  calculate  ^  in  the  same  way  of 
^  and  because  we  do  not  have  ^ .  An  alternative 
is  to  write  equations  of  moments  in  reference  system  mf. 


Mfn 


cfiw 

dmfi 


+  D 


12 


d^w 

dfi 


■2D 


d^w  \ 
dmdt  j  ’ 


d^w 

dmfi 


+ 


d^w  ^  ,  d^w 

'22^+2Z)26^j> 


(54) 

(55) 


and 


^mt 


d^w 

dm^ 


dfiw  ,  d^w\ 


(56) 


Writing  those  equations  in  a  matrix  form,  we  have: 


^m 

Mt 

,  Mffit  . 

^11  ^12  ^16 


^12  ^22  ^26 


^16  ^26  ^66 


dm^ 


d^w 

df 


d^w 
dmdt  ■ 


Thus,  from  equation  (59),  we  have: 
d^w  I  /  / 

^  ^  =  S  +  S  i2^t  +  5  jgUfmf, 
d^w 

~  ^  \2^tn  "t"  ^22^^  ^26^^^’’ 


and 


dh 


dmdt 

where 

S'  = 


S'u  S 


12  ^  )6 
^12  ^22  ^26 
^16  ^26  ^66 


=  D 


Isolating  the  unknown  variables,  we  have: 

d^w  /  /  / 

~~Q  2  ^  \2^f  ^ 


,  ,  d^w  ' 

^ 22^t  +  S —  Q^2  ~  ^12^^^ 


and 


5  26 4^^  +  S  — 


dfi\ 


dmdt 


5 


(59) 


(60) 

(61) 

(62) 

(63) 

(64) 

(65) 

(66) 
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that  can  be  written  in  the  matrix  form  as: 


-1 

‘^12 

5'i6 

f  d^w  ' 
dm- 

0 

^22 

^26 

• 

Mt 

0 

^26 

- 

.  ^mt  , 

^  _  c'  M 

af-  ^ 


d  'V  -  M 


dmdt 


(67) 


o2 

The  unknowns  M,  and  M,„,  can  be  computed  by 
solving  the  linear  system  (67). 

Finally,  the  transformation 


d^w 

'  d-w 

dx^ 

dm- 

d^w 

dy^ 

.  =T- 

d~w 

df- 

d^w 

.  dxdy 

'  dmdt  ' 

can  be  used  to  compute  second  derivatives  necessary  to 
calculate  moments  and  stresses. 


5.  Numerical  Results 

5.7.  Laminate  plate  with  clamped  edges 

Consider  a  cross-ply  laminated  graphite/epoxy  com¬ 
posite  square  plate  with  clamped  edges  under  uniformly 
distributed  load  of  intensity  q  and  with  edge  length 
a  =  1  m.  The  laminate  is  a  nine  ply  symmetrical  cross- 
ply  laminate  with  the  lay-up  [0/90/0/90/0/90/0/90/0]. 
All  plies  have  the  same  thickness.  The  total  thickness  is 
equal  to  h  -  0.001  m  and  material  properties  are:  Ei 
-  207  GPa,  Et  -  5.2  GPa,  Gu  =  3.1  GPa,  and  vu- 
0.25.  The  plate  was  discretized  using  12  quadratic  dis¬ 
continuous  boundary  elements,  as  shown  in  Figure  4. 
The  displacement  at  point  A  and  the  moment  at  point  B, 
shown  in  Figure  4,  are  compared  with  finite  element  re¬ 
sults  obtained  by  [23].  As  it  can  be  seen  in  Table  1,  the 
agreement  between  the  boundary  element  thin  plate  and 
the  finite  element  shear  deformable  plate  is  very  good  if 
we  considered  that  the  finite  element  formulation  takes 
into  account  the  effects  of  shear  deformation. 

Figure  5  shows  the  distribution  of  the  stress  ctx  at 
point  B  along  the  thickness  of  the  plate,  from  the  mid¬ 
surface  to  the  top  surface.  It  can  be  seen  that,  as  the 
stiffness  of  the  material  is  higher  in  the  direction  of  the 
fibers,  the  stress  is  higher  in  the  lamina  with  fibers  ori¬ 
ented  parallel  to  the  axis  x{a  -  0°). 


Figure  4:  Discretization  square  plate  using  3  elements  per  side. 


Table  1 :  Displacements  and  moments  in  a  square  plate  with  clamped 
edges. 


Node 

This 

work 

Reference 

[23] 

A 

wEjh^ ! 
{qa^)  X  10^ 

0.9511 

0.9341 

B 

Mxxl 

iqa^)  X  10^ 

-7.0704 

-6.6551 

5.2.  Laminated  plate  with  simply-supported  edges 

Now,  consider  the  same  laminate  with  stacking  se¬ 
quence  {+al  -aj  -\-al-al  +  al-al  +  al-al-\-a\  with 
0  <  a  <  45".  All  edges  are  simply-supported.  Material 
properties,  dimensions,  load,  and  the  mesh  used  are  the 
same  of  previous  example.  Figures  6  and  7  show  the 
effect  of  the  variation  of  9  on  the  displacement  and  re¬ 
sultant  moments,  respectively,  at  the  centre  of  the  plate 
(point  A).  They  are  compared  with  finite  element  results 
obtained  by  [23].  As  it  can  be  seen,  in  both  cases  the 
agreement  between  the  boundary  element  thin  plate  and 
the  finite  element  shear  deformable  plate  is  very  good. 
Figure  8  shows  the  stress  distribution  cr^  along  the  thick¬ 
ness  of  the  plate  at  point  A,  considering  a  -  45".  In  this 
case,  as  all  layers  present  fibers  inclined  with  respect 
to  X  direction,  variations  of  stress  between  layers  are 
smaller  than  in  previous  case. 

5.3.  Orthotropic  plate  with  clamped  edges 

The  third  example  is  a  clamped  square  plate  with  edge 
length  a  -  0.254  m  and  the  ratio  between  thicknesses 
and  edge  length  hja  -  0.05.  The  plate  is  subjected 
to  a  uniformly  distributed  static  load  q.  The  following 
material  parameters  are  used  in  the  numerical  analy¬ 
sis:  Ej  -  6.895  GPa,  Ei  -  lEj,  vu  -  0.3,  and 
Glt  -  EtI2(1  -h  vit)-  The  mesh  used  is  the  same 
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Figure  5:  Stress  distribution  along  the  thickness  for  the  cross-ply 
laminate  at  point  B  of  the  plate. 


Figure  6:  Effect  of  the  fiber  orientation  a  on  the  transversal  displace¬ 
ment  response  at  the  centre  of  the  plate. 

of  the  first  example  (Figure  4).  Figure  9  presents  mo¬ 
ments  Mx  computed  by  the  present  work  and  results 
obtained  by  [24],  using  the  meshless  Petrov-Galerkin 
method  (MLPG).  Results  are  shown  along  the  central 
line  of  the  plate,  at  y  =  a/2.  As  it  can  be  seen,  there 
is  a  perfect  agreement  between  both  results.  As  in  [24], 
bending  moments  are  normalized  by  the  central  bending 
moment  value  of  isotropic  plate  M^'‘™\a/2)  =  3064  Nm. 

6.  Conclusions 

This  paper  presented  a  boundary  integral  formula¬ 
tion  for  the  computation  of  moments  and  stresses  at 
internal  and  boundary  points  of  laminated  composite 
thin  plates.  An  integral  equation  for  the  second  dis¬ 
placement  derivative  is  developed  and  all  derivatives  of 
the  fundamental  solution  are  computed  analytically.  In 
the  proposed  approach,  in  order  to  avoid  singularities 
that  are  higher  than  hypersingular,  second  derivatives 
of  transversal  displacement  on  the  boundary  were  com¬ 


Figure  7 :  Effect  of  the  orientation  a  on  the  moment  response  at  the 
centre  of  the  plate. 


Figure  8:  Stress  distribution  (Tx  along  the  thickness  for  a  =  45^  at  the 
centre  point  of  the  plate. 

puted  using  constitutive  equations  and  shape  function 
derivatives.  The  obtained  results  are  in  good  agreement 
when  compared  with  results  available  in  literature. 
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